TransitionMatrixPackage

Yet To Be Implemented ...

For Version 1.0 Beta

Error in MFPT

Straight Forward, easy to implement

Hummer Method to get Rate Matrices

Read, how it works and then try to implement it
- Easy concept, but not easy to implement with simple Mathematica Functions
- Got solutions for a 4x4 system...
- Needs a more simple method to get solution. Maybe possible to adapt the Boxer Method to result in same solution.

Franks Method to sample Matrices from Dirichlet Distribution with Detailed Balance Constraint, Positivity Constraint and/or given Equilibrium Distribution

Should be easy, but not sure yet, what is the easiest way to add it to the library, maybe also as Distribution, without Mean / Variance, just DrawTransitionMatrix / DrawTransitionMatrixSeries implemented, Options may be given to set initial matrix

Visualization

Use some things from other former notebooks. Maybe make a suitable colorrange. Makes no sense to incluse 0.5 in the colorrange, since this is not precise enough. Use Mesh -> 0.5 instead.

CleanUp

Error Messages

Documentation

Nothing Yet, keep it short and easy to understand. Give lots of examples

For Further Versions

Error in TPT !!! for Frank

Still need to write down the calculations

Helpers for Sequences of 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&]]

Options

Options[TransitionMatrix]={Normalized->True,ShowSize->7,Lagtime->1};

Options[RateMatrix]={Normalized->True,ShowSize->7};

Options[CountMatrix]={Distribution->DirichletDistribution,Prior->1,TransitionMatrix->None,ShowSize->7};

Options[TimeSeries]={States->1,Lagtime->1,ShowSize->11};

Options[CrispClusterMatrix]={Complete->False,States->0};

Options[KMeansClustering]={Iterations->20};

Options[ClusterMatrix]={};

Options[ToRateMatrix]={Method->MatrixLog,Lagtime->Inherit};

Options[DirichletWithDBDistribution]={Iterations->1000,MaximumResidue->0.001};

Message Definitions

MaxLikelihoodMatrix::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.";

Matrixoperationen

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

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

CountMatrices

CountMatrix/:MakeBoxes[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],StandardForm]:=If[OptionValue[CountMatrix,ShowSize]>Length[m],InterpretationBox[#,CountMatrix[m,opts]]&[RowBox[{SubscriptBox["ℂ",Length[m]],ToBoxes[MatrixForm[m]]}]],InterpretationBox[#,CountMatrix[m,opts]]&[SubscriptBox["ℂ",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/:MaxLikelihoodMatrix[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],optsA___]:=MaxLikelihoodMatrix[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]]

CountMatrix/:CountMatrix[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],optsCM:OptionsPattern[TransitionMatrix]]:=Module[{cM,options,nStates,nLagtime},
options=UnionRules[opts,States->Max[v]];
nStates=Max[OptionValue[TimeSeries,States],Max[v]];
nLagtime=OptionValue[TimeSeries,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]]]

RateMatrix Definitions

RateMatrix/:MakeBoxes[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]],StandardForm]:=If[OptionValue[RateMatrix,ShowSize]>Length[m],InterpretationBox[#,RateMatrix[m,opts]]&[RowBox[{SubscriptBox["ℝ",Length[m]],ToBoxes[MatrixForm[m]]}]],InterpretationBox[#,RateMatrix[m,opts]]&[SubscriptBox["ℝ",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/: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]]]:=OptionValue[RateMatrix,Normalized]

RateMatrix/:Normalize[RateMatrix[m_?MatrixQ,opts:OptionsPattern[RateMatrix]]]:=RateMatrix[Map[#-Apply[Plus,#]&,m],ReplaceRules[Normalized->True,opts]]

RateMatrixFromMatrix[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 Definitions

TransitionMatrix/:MakeBoxes[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],StandardForm]:=If[OptionValue[TransitionMatrix,ShowSize]>Length[m],InterpretationBox[#,TransitionMatrix[m,opts]]&[RowBox[{SubscriptBox["𝕋",Length[m]],ToBoxes[MatrixForm[m]]}]],InterpretationBox[#,TransitionMatrix[m,opts]]&[SubscriptBox["𝕋",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

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/: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]]

TransitionMatrixFromMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]:=TransitionMatrix[Map[#/Apply[Plus,#]&,m],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,Length[states]];
Inverse[cMEV[[All,states]]].cMEV
]

TransitionMatrix/:Equilibrium[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]]]:=#/Apply[Plus,#]&[Eigenvectors[Transpose[m],1][[1]]]

TransitionMatrix/:ToRateMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],optsTRM:OptionsPattern[ToRateMatrix]]:=Module[{rM},
If[OptionValue[ToRateMatrix,optsTRM,Method]===MatrixLog,
{eVals,eVecs}=Eigensystem[m];
nLagtime=If[OptionValue[ToRateMatrix,optsTRM,Method]===Inherit,
OptionValue[TransitionMatrix,Lagtime],
OptionValue[ToRateMatrix,optsTRM,Lagtime]
];
rM=RateMatrix[Transpose[eVecs].DiagonalMatrix[Log[eVals]/nLagtime].Inverse[Transpose[eVecs]],Lagtime->nLagtime];
];
If[OptionValue[ToRateMatrix,List[optsTRM],Method]===Hummer,
Message[ToRateMatrix::MethodNotImplemented];
rM=Null;
];
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]]]]]]]
]

TransitionMatrix/:MeanFirstPassageTime[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_]:=Module[{cM,cMEV},
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]}]
]

TransitionMatrix/:EquilibriumSensitivityMatrix[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/:CommittorSensitivityVector[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[[Complement[Range[size],states],Complement[Range[size],states]]]-IdentityMatrix[size-Length[states]]];
cA=Transpose[cS].cM.cS;
Table[Apply[Plus,cA[[ii]]],{ii,1,size}]
]

TransitionMatrix/:ChooseTransition[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],state_,num___:1]:=Module[{probabilities,probSum,endState,p,size,list,startState},
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/: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]]]

Distribution Handling

DirichletDistribution/:MeanMatrix[DirichletDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]]:=TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m+OptionValue[CountMatrix,Prior])]]

DirichletDistribution/:MaxLikelihoodMatrix[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}]]]

FixedPointDistribution/:MeanMatrix[FixedPointDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]]:=TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m)]]

FixedPointDistribution/:MaxLikelihoodMatrix[FixedPointDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]]:=TransitionMatrixFromMatrix[Map[#/Apply[Plus,#]&,(m)]]

FixedPointDistribution/:VarianceMatrix[FixedPointDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]]:=0*m

FixedPointDistribution/:DrawTransitionMatrix[FixedPointDistribution[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/:MaxLikelihoodMatrix[DirichletWithDBDistribution[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]],optsML:OptionsPattern[DirichletWithDBDistribution]]:=
Module[{size,sumVec,symmMat,result,nIterations,sumResVec,nResidue},
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],MaximumResidue];
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]]
]

Trajectories / Timeseries

TimeSeries/:MakeBoxes[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],StandardForm]:=
If[Length[v]<OptionValue[TimeSeries,ShowSize],InterpretationBox[#,TimeSeries[v,opts]]&[RowBox[{SubscriptBox["𝕊",RowBox[{Max[OptionValue[TimeSeries,States],Max[v]]}]],"(",ToBoxes[TableForm[v,TableDirections->Row]],")"}]],InterpretationBox[#,TimeSeries[v,opts]]&[RowBox[{SubscriptBox["𝕊",RowBox[{Max[OptionValue[TimeSeries,States],Max[v]]}]],"⟨",Length[v],"⟩"}]]]

TimeSeriesQ[x_]:=If[Head[x]==TimeSeries&&VectorQ[x[[1]]],True,False];

TimeSeries/:NormalForm[TimeSeries[v_?MatrixQ,opts:OptionsPattern[TimeSeries]]]:={v,opts}

TimeSeries/:Normal[TimeSeries[v_?MatrixQ,opts:OptionsPattern[TimeSeries]]]:=v

TimeSeries/:MatrixForm[TimeSeries[v_?MatrixQ,opts:OptionsPattern[TimeSeries]]]:=MatrixForm[v]

TimeSeries/:Length[TimeSeries[v_?MatrixQ],opts:OptionsPattern[TimeSeries]]:=Length[v]

Trajectory/:MakeBoxes[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],StandardForm]:=
InterpretationBox[#,Trajectory[m,opts]]&[RowBox[{SubscriptBox["𝕄",Dimensions[m][[2]]],"⟨",Dimensions[m][[1]],"⟩"}]]

Trajectory/:ToTimeSeries[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],centers_]:=Module[{},
nearest=Nearest[centers->Automatic];
TimeSeries[Map[nearest[#,1][[1]]&,m],States->Length[centers]]
]

Trajectory/:KMeansClustering[Trajectory[m_?MatrixQ,opts:OptionsPattern[Trajectory]],nCenters_,optsKM:OptionsPattern[KMeansClustering]]:=Module[{centerList,nearest,clusters,nClusters,nIterations,oCL,residuum},
centerList=If[
NumberQ[nCenters],
m[[RandomInteger[{1,Length[m]},nCenters]]],
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 & Clustering

ClusterMatrix/:MakeBoxes[ClusterMatrix[m_?MatrixQ,opts:OptionsPattern[ClusterMatrix]],StandardForm]:=InterpretationBox[#,ClusterMatrix[m,opts]]&[SubscriptBox["𝕃",RowBox[{Dimensions[m][[2]],"->",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"]

CountMatrix/:Cluster[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]]:=CountMatrix[c.m.Transpose[c],ReplaceRules[States->Dimensions[c][[1]],opts]]

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},
cM=Select[Table[{nn,Apply[Plus,m[[nn]]]},{nn,1,Length[m]}],#[[2]]>0&][[All,1]];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]}]]
]

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]}]]
]

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]]
];

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]]]}]
]

TimeSeries/:Cluster[TimeSeries[v_?VectorQ,opts:OptionsPattern[TimeSeries]],ClusterMatrix[c_?MatrixQ,optsCM:OptionsPattern[ClusterMatrix]]]:=Module[{size,cM,rules},
size=Length[m];
rules=ToCrispRules[ClusterMatrix[c,optsCM]];
v/.rules
];

CountMatrix/:HasEmptyStatesQ[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]:=If[Apply[Times,Map[Apply[Plus,#]&,m]]==0,True,False]

PCCA Clustering

TransitionMatrix/:PCCASoftClusterMatrix[TransitionMatrix[m_?MatrixQ,opts:OptionsPattern[TransitionMatrix]],level_]:=Module[{eqDist,size,cM,eVecs,distTable,maxDist,cornerList,eVecsNorm,list2},
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]];
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},
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]];
ClusterMatrix[Transpose[Crispify[Transpose[Chop[Join[Transpose[list2],{Map[1-Apply[Plus,#]&,list2]}]]]]]]
];

Flows & Transition Path Theory

TransitionMatrix/:TransitionFlowMatrix[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[[jj]]*(1-cM[[ii]])*eqDist[[ii]]*m[[ii,jj]],{ii,1,size},{jj,1,size}]
];

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}]
];

Plotfunctions

PositionDensityPlot[positionList_,valueList_,opts:OptionsPattern[DensityPlot]]:=Module[{arrayToPlot},
arrayToPlot=If[NumberQ[positionList],
Partition[valueList,positionList],
Transpose[Join[Transpose[positionList],{valueList}]]
];
ListDensityPlot[arrayToPlot,BaseStyle->{FontFamily->"GillSans",FontSize->12},opts]
]

PositionStreamDensityPlot[positionList_,valueList_,opts:OptionsPattern[Join[{FieldStyle->{Gradient}},Options[StreamPlot]]]]:=Module[{arrayToPlot,interpolation,func},
arrayToPlot=If[NumberQ[positionList],
Transpose[Join[Transpose[pos30=Table[#&[{Mod[ii,positionList]+1,Floor[ii/positionList]+1}],{ii,0,positionList^2-1}]],{valueList}]],
Transpose[Join[Transpose[positionList],{valueList}]]
];
interpolation=Interpolation[arrayToPlot];
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[FieldStyle,opts]],BaseStyle->{FontFamily->"GillSans",FontSize->12}]
]

PositionStreamPlot[positionList_,valueList_,opts:OptionsPattern[Join[{FieldStyle->{Gradient}},Options[StreamPlot]]]]:=Module[{arrayToPlot,interpolation,func},
arrayToPlot=If[NumberQ[positionList],
Transpose[Join[Transpose[pos30=Table[#&[{Mod[ii,positionList]+1,Floor[ii/positionList]+1}],{ii,0,positionList^2-1}]],{valueList}]],
Transpose[Join[Transpose[positionList],{valueList}]]
];
interpolation=Interpolation[arrayToPlot];
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[FieldStyle,opts]],BaseStyle->{FontFamily->"GillSans",FontSize->12}]
]

I/O Function

ReadSparseMatrix[fileName_,opts:OptionsPattern[Import]]:=Module[
{readImport},
readImport=Import[fileName,"Table",opts];
Normal[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]
]

TestMatrices

m3=N[{{0.9,0.00,0.1},{0.0,0.9,0.1},{0.08,0.07,0.85}}];

c2=CountMatrix[{{2,1},{1,3}}];

t2=TransitionMatrixFromMatrix[{{2,1},{1,3}}];

c4=CountMatrix[(({{20, 6, 0, 0}, {3, 50, 7, 0}, {0, 4, 80, 9}, {0, 0, 2, 30}}))];

c42=CountMatrix[(({{20, 6, 2, 0}, {3, 50, 7, 0}, {0, 4, 80, 9}, {0, 0, 2, 30}}))];

c6=CountMatrix[(({{79, 3, 4, 6, 7, 7}, {3, 93, 4, 3, 7, 5}, {6, 5, 106, 1, 3, 1}, {1, 4, 5, 88, 1, 3}, {4, 7, 5, 8, 78, 8}, {9, 7, 6, 1, 3, 86}}))];

c6=CountMatrix[TimeSeries[RandomInteger[{1,6},1000000],Lagtime->1]];

c62=CountMatrix[Normal[c6]+Transpose[Normal[c6]]];

readImport=Import["/Users/jan-hendrikprinz/Desktop/Committor Paper/data/02.t","Table"];
tM=TransitionMatrix[Normal[SparseArray[Map[1+#[[{1,2}]]->#[[3]]&,readImport]]]];
pos30=Table[#&[{Mod[ii,30]+1,Floor[ii/30]+1}+0*{Random[],Random[]}],{ii,0,899}];

t6=Normalize[CountMatrix[TimeSeries[RandomInteger[{1,6},10000000],States->6,Lagtime->2]]];

Tests

Testing

c3=CountMatrix[(({{20, 3, 4}, {1, 30, 2}, {3, 4, 25}})),Prior->1]

Options[c3]

c3DB=CountMatrix[c3,Distribution-> DirichletWithDBDistribution]

c3DB

N[MaxLikelihoodMatrix[c3]]

MaxLikelihoodMatrix[c3DB]

c42DB=CountMatrix[c4,Distribution-> DirichletWithDBDistribution]

N[Normalize[c42DB]]

MaxLikelihoodMatrix[c42DB,Iterations->100000,MaximumResidue->0.0000001]

mat=UpperTriangularize[Normal[c4]]//MatrixForm

mat=Table[t[ii,jj],{ii,1,4},{jj,1,4}]

mat2=LUDecomposition[mat];

mat2=LUDecomposition[mat];
UDmat=UpperTriangularize[mat2]
LDmat=LowerTriangularize[mat2]

MatrixExp[A.D]==MatrixExp[A]+MatrixExp[D]

MaxLikelihoodMatrix[c42DB]

m2=Table[Random[],{ii,1,3},{jj,1,3}];
v2=Table[Random[],{jj,1,3}];

Eigenvalues[m2+DiagonalMatrix[{0,0,0}]]
Eigenvalues[m2]

MatrixExp[m2.DiagonalMatrix[v2]]-MatrixExp[m2]

Transpose[DiagonalMatrix[{p[1],p[2],p[3],p[4]}].UpperTriangularize[mat].Inverse[DiagonalMatrix[{p[1],p[2],p[3],p[4]}]]]+UpperTriangularize[mat]

MatrixForm[mat4]

Transpose[DiagonalMatrix[{p[1],p[2],p[3],p[4]}].UpperTriangularize[mat].Inverse[DiagonalMatrix[{p[1],p[2],p[3],p[4]}]]]+UpperTriangularize[mat]

CountMatrix/:Hummer[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]:=Module[{size,cons,expression,vars,mat},
size=Length[m];
mat=MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]];
expression=Sum[m[[ii,jj]]*Log[mat[[ii,jj]]],{ii,1,size},{jj,1,size}];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
vars=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];
NMinimize[{-Re[expression],cons},vars,MaxIterations->1000,PrecisionGoal->20]
]

mat4=Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]/.size->4;

mat3=Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]/.size->3;

mat3

Eigenvalues[mat4]//FullSimplify

func=Apply[Function,{vars,
If[NumberQ[p[1]],MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]],IdentityMatrix[size]]}];
func2=Apply[Function,{vars,Sum[m[[ii,jj]]*Apply[func,vars][[ii,jj]],{ii,1,size},{jj,1,size}]}];

mat4log=Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]/.size->4;

CountMatrix/:Hummer[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]:=Module[{size,cons,expression,vars,mat,varsSet},
size=Length[m];
vars=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[p[kk],{kk,1,size-1}]];
func=Apply[Function,{vars,
If[NumberQ[p[1]],MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]],IdentityMatrix[size]]}];
func2=Apply[Compile,{vars,Apply[Function,{vars,Sum[m[[ii,jj]]*Re[Apply[func,vars][[ii,jj]]],{ii,1,size},{jj,1,size}]}]}];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
varsSet=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];
NMinimize[{Apply[func2,vars],cons},varsSet]
]

CountMatrix/:Hummer[CountMatrix[m_?MatrixQ,opts:OptionsPattern[CountMatrix]]]:=Module[{},
size=Length[m];
vars=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[p[kk],{kk,1,size-1}]];
func=Apply[Function,{vars,
If[NumberQ[p[1]],MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]],IdentityMatrix[size]]}];
func2=Apply[Compile,{vars,Apply[Function,{vars,Sum[m[[ii,jj]]*Re[Apply[func,vars][[ii,jj]]],{ii,1,size},{jj,1,size}]}]}];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
varsSet=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];
NMinimize[{Apply[func2,vars],cons},varsSet]
]

func2

vars

func2

g[lst:{t_,r_}/;And@@NumericQ/@lst]:=
(x1-1)^2+(x2-1)^2-x3^2+(x4+1)^2+x5^2

func2[{0.2,0.2,0.2,0.2,0.2}]

g[{0.1,0.2}]

m={{1,3,2},{3,2,3},{4,1,1}};
size=Length[m];
vars=Sequence[Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[p[kk],{kk,1,size-1}]]];
mat=Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}];
func2:=Function[{vars},Sum[m[[ii,jj]]*Re[Log[MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]][[ii,jj]]]],{ii,1,size},{jj,1,size}]];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
varsSet=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];

func2[{0.2,0.2,0.2,0.2,0.2}]

Remove[func2]

NMinimize[{Apply[func2,vars],cons},varsSet]

?func2

fmat

func2

ReleaseHold

Hummer[CountMatrix[{{1,3,3},{2,3,1},{1,2,4}}]]

{evals,evecs}=Eigensystem[N[Normal[Normalize[c4]]]];

matc4=Chop[Transpose[evecs].DiagonalMatrix[evals].Transpose[Inverse[evecs]]]

matc4log=Chop[Transpose[evecs].DiagonalMatrix[Log[evals]].Transpose[Inverse[evecs]]]

vars=Join[Flatten[Table[t[ii,jj],{ii,1,6-1},{jj,ii+1,6}]],Table[p[kk],{kk,1,6-1}]];
func=If[NumberQ[p[1]],MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==6,1-Sum[p[ll],{ll,1,6-1}],p[ii]],t[ii,kk]],{kk,1,6}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==6,1-Sum[p[kk],{kk,1,6-1}],p[ii]],1],{ii,1,6},{jj,1,6}]],IdentityMatrix[6]];

Module[{m,size,cons,expression,vars,mat,matExp,sol},
m=Normal[c3];
size=Length[m];
mat=MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]];
matExp=MatrixExp[mat];
expression=Sum[m[[ii,jj]]*Log[matExp[[ii,jj]]],{ii,1,size},{jj,1,size}];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
vars=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];
sol=NMinimize[{Abs[expression],cons},vars];
mat/.sol[[2]]
]

Module[{size,cons,vars},
mat=N[Normal[c4]];
size=Length[mat];
cons=Apply[And,Flatten[Join[Table[p[kk]≥0,{kk,1,size-1}],{Sum[p[kk],{kk,1,size-1}]≤1},Table[t[ii,jj]≥0,{ii,1,size-1},{jj,ii+1,size}]]]];
vars=Join[Flatten[Table[t[ii,jj],{ii,1,size-1},{jj,ii+1,size}]],Table[{p[kk],0.8/size,1.2/size},{kk,1,size-1}]];
NMinimize[{-Abs[Sum[MatrixExp[Table[If[ii==jj,t[ii,ii]-Sum[If[ii>kk,t[kk,ii]*p[kk]/If[ii==size,1-Sum[p[ll],{ll,1,size-1}],p[ii]],t[ii,kk]],{kk,1,size}],If[ii<jj,t[ii,jj],t[jj,ii]]]*If[ii>jj,p[jj]/If[ii==size,1-Sum[p[kk],{kk,1,size-1}],p[ii]],1],{ii,1,size},{jj,1,size}]][[ii,jj]]*Log[mat[[ii,jj]]],{ii,1,size},{jj,1,size}]],cons},vars]
]

MeanMatrix[c42]

c42

ImpliedTimeScales[t4]

Eigenvectors[Transpose[Normal[ToRateMatrix[t4]]]][[-1]]

Equilibrium[t4]

Equilibrium[ToRateMatrix[t4]]

Chop[Normal[ToTransitionMatrix[ToRateMatrix[t4],3]]-Normal[t4]]

RemoveClusterMatrix[{1,2,3},10]

Normal[RemoveEmptyStatesClusterMatrix[c6]]

Cluster[c6,CrispClusterMatrix[{{1,2,3},{4,5},{5,6},{}}]]

Cluster[c6,{{0,.2,1,1,.2,0},{1,.8,0,0,.8,1}}]

pccaMembership=PCCASoftClusterMatrix[tM,3];

Map[Apply[Plus,#]&,Transpose[pccaMembership]]

tM3=Cluster[tM,pccaMembership];

N[TransitionMatrix[MeanMatrix[c6],States->10]]

tM3

reactivity=ReactivityVector[tM,280,640];

plot1=PositionStreamPlot[pos30,Log[reactivity+0.00000000001],FieldStyle->{Rotation},StreamStyle->Directive[{Black}],StreamScale->None,StreamPoints->Fine]

plot2=PositionDensityPlot[30,Log[reactivity],InterpolationOrder->1,ColorFunction-> "TemperatureMap",RegionFunction->((Norm[{#1,#2}-{10,10}]>Sqrt[2])&&(Norm[{#1,#2}-{10,22}]>Sqrt[2])&),PlotPoints->300,Mesh->1,MeshFunctions->{#3&}]

Show[plot2,plot1]

Committor[N[MeanMatrix[c6]],{2}]

tfM=TransitionFlowMatrix[tM,280,640];

Apply[Plus,Flatten[tfM]]

Apply[Plus,Flatten[tfM]]

Apply[Plus,tfM[[All,640]]]

Apply[Plus,tfM[[280]]]

symmTfM=Map[Max[#,0]&,tfM-Transpose[tfM],{2}];

eq=Equilibrium[tM];

Nearest[pos30->Automatic,{10,22}]

comm1=Committor[tM,{280,640}];

plot1=PositionStreamPlot[30,comm1[[1]],FieldStyle->{Gradient},StreamStyle->Black,StreamPoints->Fine];
plot2=PositionDensityPlot[30,comm1[[1]],InterpolationOrder->3,ColorFunction->"TemperatureMap",Mesh->{{0.01,0.5,0.999}},MeshFunctions->{#3&},MeshStyle->{Thick}];
Show[{plot2,plot1}]

plot1=PositionStreamPlot[pos30,pccaMembership[[1]],FieldStyle->{Gradient},StreamStyle->Directive[{Black}]];
plot2=PositionDensityPlot[pos30,pccaMembership[[1]],InterpolationOrder->2,ColorFunction->"TemperatureMap"];
Show[{plot2,plot1}]

PositionDensityPlot[30,Log[eq],InterpolationOrder->2,ColorFunction->"TemperatureMap"]

plotLevel=1/2;
Show[Table[
PositionDensityPlot[pos30,pccaMembership[[nn]],Mesh->{{plotLevel}},MeshFunctions->{Function[{x,y,z},z]},MeshStyle->Directive[Thick,Black],ColorFunction->({Opacity[#],ColorData[35,"ColorList"][[nn*3]]}&),RegionFunction->Function[{x,y,z},z>plotLevel*0.999],PlotRange->{0,1}],
{nn,1,3}]]

Manipulate[Show[Table[
pos30=Table[{Floor[ii/30]+1,Mod[ii,30]+1},{ii,0,899}];
arrayToPlot=Transpose[Join[Transpose[pos30],Transpose[Map[#[[{nn}]]&,Transpose[list]]]]];
ListDensityPlot[arrayToPlot,Mesh->{{plotLevel}},MeshFunctions->{Function[{x,y,z},z]},MeshStyle->Directive[Thick,Black],BaseStyle->{FontFamily->"GillSans",FontSize->12},FrameTicks->All,GridLines->Automatic,ColorFunction->(RGBColor[If[nn==1,1,0],If[nn==2,1,0],If[nn==3,1,0],0.1#]&),RegionFunction->Function[{x,y,z},z>plotLevel*0.95]],
{nn,1,3}]],{plotLevel,0,1}]

ArrayPlot[Partition[Transpose[list],30],ColorFunction->RGBColor,Frame->True,FrameTicks->All,DataReversed->True,BaseStyle->{FontFamily->"GillSans",FontSize->12}]

PositionDensityPlot[30,eq]

statePlot=Graphics[MapIndexed[{GrayLevel[0],Disk[#,.707],GrayLevel[1],Style[Text[FromCharacterCode[#2[[1]]+64],#+{0,0}],FontFamily->"GillSans"]}&,Map[(.5+GetPositionFromStateC[#])&,cclist]]];

MatrixPlot[Partition[eq,30],BaseStyle->{FontFamily->"GillSans",FontSize->12}]

Normal[c6]
N[Normalize[c62]]
N[Symmetrize[Normalize[c6]]]
N[Symmetrize[Normalize[c62]]]

s=ChooseTransition[t6,1,100000]

Normalize[CountMatrix[s]]//N

Normal[Normalize[CountMatrix[s]]]//N//MatrixForm
Normal[t6]//N//MatrixForm
(Normal[Normalize[CountMatrix[s]]]-Normal[t6])//N//MatrixForm

CountMatrix[TimeSeries[RandomInteger[{1,6},100],States->6,Lagtime->50]]

tab1=Table[GetInstance[c4]-Normal[MeanMatrix[c4]],{100000}];

Variance[tab1]-N[VarianceMatrix[c4]]

Normal[MeanMatrix[c2]]

VarianceMatrix[c2]

N[CommittorSensitivityVector[MaxLikelihoodMatrix[c4],{1,4}]]

Apply[Plus,Flatten[N[Abs[EquilibriumSensitivityMatrix[MaxLikelihoodMatrix[c4]]]]*VarianceMatrix[c4],{1,2}]]

cN=Normalize[c2]

NormalForm[TransitionMatrix[MeanMatrix[c]]]

Normal[Normal[Normalize[TransitionMatrix[SparseArray[{{1,1}->1,{2,2}->2}]]]]]

Normal[Normalize[c]]