Loop Eigenvalue Elasticity Analysis (LEEA) Utility -- Instructions

Purpose
This Mathematica notebook contains all routines to perform the Loop Eigenvalue Elasticity Analysis (LEEA) as described in Kampmann and Oliva (2006).

Requirements
The LEEA Utility notebook requires that the FeedbackLoops.m package is loadable into Mathematica.  The feedbackLoops.m package is available from http://iops.tamu.edu/faculty/roliva/research/sd/leea/toolset.html.  See section A.8.2 of the Mathematica Book for directory locations for the file to be loadable to Mathematica.
 
Input files
	1. A Mathematica notebook file (*.nb) that contains the structure of the model that will be analyzed.  
		An utility to translate Vensim models into this format is available at http://nico.tamu.edu/mdl2nb/mdl2nb.html.  Be sure to address the limitations of the mdl2nb utility before using the LEEA Utility.
	2. A tab delimited file with the model simulation output.  
		Use Vensim Model>Dataset Export ... command (with options "tab" and "time running down") to translate the Vensim data file (*.vdf) into a tab delimited file (*.tab).  The LEEA Utility will calculate the LEEA for every Time Interval (dt) stored in the *.tab file.  Use Vensim Saveper option to control the frequency of calculations.
	
Output files
	1. A tab delimited file containing the Loop Elasticities and Influences. 
		The file contains, in a table format, the evolution (for every Saveper included in the *.tab input file) of the real and imaginary parts of the eigenvalues and the loop gains of every loop in the model's Shortest Independent Loop Set (SILS). Additionally, the file contains the real, imaginary and absolute values of elasticity and influence of every loop in the model's SILS to every eigenvalue for every Saveper included in the *.tab input file.  The file can be imported into Excel to report and evaluate the results of the LEEA.
	2. This notebook itself will take the values from the computations with the specified model and as such it can be saved to retain partial results (e.g., adjacency matrix, Shorter independent loop set, model Jacobian, etc.).
	
Instructions
1. Before getting started it is useful to turn off the message General::spell (Off[General::spell];) to avoid harmless, but annoying, warnings that might be triggered by Vensim variable names. 
2. Open (File>Open ...) the notebook file (*.nb) that contains the structure of the model that will be analyzed and evaluate it (Kernel>Evaluation>Evaluate Notebook).
3. Bring to the forefront of Mathematica the LEEA.nb notebook (this notebook) and evaluate it ((Kernel>Evaluation>Evaluate Notebook).   
	The utility will first ask for the location of the *.tab file with the model simulation output.  Once the file is located and selected, hit <<Open>> and the utility will proceed to do the computations.  The utility will automatically save two tab delimited files, containing Link and Loop Elasticities in the current default directory (use Directory[] to identify current default directory and locate the files).  
4. Save this notebook with a different name (File>Save as ...) to retain partial results of the LEEA analysis (e.g., adjacency matrix, shorter independent loop set, model Jacobian, etc.).

Limitations
This implementation of LEEA works only with continuous (differentiable) table functions that are expressed as functions among the variable list.  Please ensure that the Vensim model that is imported into Mathematica adheres to this requirement. 

History
v0.1	02/2006		Core functionality with explicit model and data inputs.
v0.2	02/2006		Eliminated the link elasticity report and included the influence computations.

(c) Christian E. Kampmann and Rogelio Oliva, 2006. All Rights Reserved.

Import simulation data

Obtaining simulation data file

Needs["GUIKit`"];
\(\(FileName = 
        GUIResolve[
Script[
Widget["FileDialog", 
              Name -> "openFileDialog"]; 
SetPropertyValue[{"openFileDialog", "multiSelectionEnabled"}, \ 
              False]; 
SetPropertyValue[{"openFileDialog", "fileSelectionMode"}, \ 
PropertyValue[{"openFileDialog", "FILES_ONLY"}]]; 
returnValue = 
              InvokeMethod[{"openFileDialog", "showOpenDialog"}, \ 
                Null]; 
If[
              returnValue\  === 
                PropertyValue[{"openFileDialog", "APPROVE_OPTION"}], 
              PropertyValue[{PropertyValue[{"openFileDialog", "selectedFile"}], \ "path"}], 
Null
]
]
];\)\)
imp = Import[FileName, "TSV"];
imp[[1]] = 
      StringReplace[imp[[1]], " " -> ""]

Data preparation

For[i = 3,  
        i ≤ Dimensions[
              imp][[1]], i++, 
For[j = 1,  
            j ≤ Dimensions[
                  imp][[2]],  j++, 
If[
                imp[[i, j]] == "", 
imp[[i, 
                      j]] = 
                  imp[[2, 
                      j]]
];];];
SubNum[n_] := 
      MapThread[
        Rule[ToExpression[#1], ToExpression[#2]] &, {StringReplace[
            imp[[1]], " " -> ""], imp[[n]]}]

Model as graph

Needs["FeedbackLoops`"];

Nodes

nodelist=Join[xv,yv];

nodelist

Adjacency matrix and strong components

Mm=FromAdjacencyLists[
	Block[
		{t=Join[dsv,aux]},
		Map[(Map[First,Position[t,#]])&,nodelist]
	]
];

sc = StrongComponents[Mm]

Edges

el=ToEdgeList[Mm]

edgelist = Map[nodelist[[#]] &, el, {2}]

Directed cycles and cycle matrix

cl = DirectedCycles[Mm];

clg=GeodeticCycles[Mm];

cm = CycleMatrix[cl, el];

cmg = CycleMatrix[clg, el];

clText = Map[Function[nodelist[[#1]]], clg, {2}];

Statistics and checks

noesc = 
        Length[Cases[
            Map[First[First[Position[sc, #]]] &, el, {2}], {x_, 
              x_}]];

Print["No. of edges:            ",Length[el]];
Print["No. of strong comp.:     ",Length[sc]];
Print["No. of edges in str. c.: ",noesc];
Print["No. of nodes:            ",Length[Mm]];
Print["No. of loops:            ",Length[cl]];
Print["No. of geodectic loops:  ",Length[clg]];
Print["No. of indep. fb. loops: ",noesc-Length[Mm]+Length[sc]];
Print["Check: rank of Cm:       ",Rank[cm]];
Print["Check: rank of Cmg:      ",Rank[cmg]];

A shorter independent loop set

Shorter independent loop set

clr = SILS[clg]

Cycle matrix (reduced)

cmr = CycleMatrix[clr, el];

ILS list

loopsetText = 
      Table[{"Loop " <> ToString[i], 
          Map[nodelist[[#]] &, clr[[i]]]}, {i, 
          Length[clr]}]

Statistics and checks

Print["*** Check that it is OK ***"];
Print["Every loop in original cycle? ", 
        Union[Map[MemberQ[cl, #] &, clr]]];
Print["Dimensions of cmr             ", 
        Dimensions[cmr]];
Print["Rank of cmr                   ", Rank[cmr]];
Print["Rank of cm                    ", Rank[cm]];

Expressing loop and edge gains

List edge gains g[i,j] and loop gains g[i]

eg = Map[Apply[gain, #] &, el]

      
lgr = Table[gain[i], {i, Length[cmr]}]

Edge gains as functions of variables

egv=Map[
	Apply[
		Function[
			{x1,x2},
			If[x2<=n,
				If[x1<=n,
					D[dsv[[x2]],xv[[x1]]],
					D[dsv[[x2]],yv[[x1-n]]]
				],
				If[x1<=n,
					D[aux[[x2-n]],xv[[x1]]],
					D[aux[[x2-n]],yv[[x1-n]]]
				]
			]
		],
		#
	]&,
	el
] (* Edge gains, as functions of variables *)

subEgToEgv=Table[
        eg[[i]] -> egv[[i]], {i, Length[eg]}]

Loop gains

lg = cmr . eg /. Plus -> Times;

Jacobian

Am=Table[If[MemberQ[edgelist,{xv[[j]],xv[[i]]}],gain[j,i],0],{i,n},{j,n}];
Bm=Table[If[MemberQ[edgelist,{yv[[j]],xv[[i]]}],gain[j+n,i],0],{i,n},{j,p}];
Cm=Table[If[MemberQ[edgelist,{xv[[j]],yv[[i]]}],gain[j,i+n],0],{i,p},{j,n}];
Dm=Table[If[MemberQ[edgelist,{yv[[j]],yv[[i]]}],gain[j+n,i+n],0],{i,p},{j,p}];
Um=Inverse[IdentityMatrix[p]-Dm];
Jm=Am+Bm.Um.Cm//Expand;

Link/loop elasticities

Matrix method

outLoopEla = {{"t", "i", "Re[z]", "Im[z]", "Loop nr\
", "Loop", "Nodes", "gain", "Re[e]", "Im[e]", "Abs[e]\
",  "Re[I]",  "Im[I]",  "Abs[I]"}};

Block[{x, y, i, j}, 
For[loopnodes = {}; i = 1, 
          i ≤ Length[loopsetText], i++, 
x = 
            Map[ToString, loopsetText[[i, 2]]]; 
y = 
            x[[1]]; 
For[j = 2, 
            j ≤ Length[x], j++, 
            y = y <> ">" <> 
                ToString[x[[j]]]]; 
AppendTo[loopnodes, 
            y]
]
];

\(\(Block[
{i, j, k, s, t, Ax, Bx, Cx, Dx, lx, rx, An, 
          Bn, Cn, Dn, Un, Jn, cmrinv, 
ev, es, Rn, Ln, 
          edgeGain, loopGain, linkEla, loopEla, z, \ 
          linkElaMa}, 
cmrinv = 
          N[PseudoInverse[Transpose[cmr]]]; 
linkElaMa = 
          MapThread[\((#1\ #2/
                    z)\) &, 
{Map[
\({Ax, 
                        Bx, Cx, 
                        Dx}[\([#[\([1]\)]]\)]\)[#[\([2]\)], #[\([3]\)]]\  &, 
Map[Flatten[Position[{Am, Bm, Cm, Dm}, #]] &, 
                  eg]
], 
eg
}
]; 
Ax[i_, j_] := 
          lx[\([i]\)] rx[\([j]\)] An[\([i, j]\)]; 
Bx[i_, 
            j_] := lx[\([i]\)] gx[\([j]\)]; 
Cx[i_, j_] := 
          fx[\([i]\)] rx[\([j]\)]; 
Dx[i_, j_] := 
          fx[\([i]\)] gx[\([j]\)]; 
For[i = 2, 
          i ≤ Length[imp], \(i++\), 
s = 
            subEgToEgv /. SubNum[i]; 
t = 
            Time /. SubNum[i]; 
Print[
            t]; 
An = Am /. s; 
Bn = 
            Bm /. s; 
Cn = Cm /. s; 
Dn = 
            Dm /. s; 
Un = Um /. s; 
Jn = 
            Jm /. s; 
es = 
            Eigensystem[Jn]; 
ev = 
            es[\([1]\)]; 
Rn = 
            es[\([2]\)]; 
Ln = 
            Inverse[Transpose[Rn]] // Chop; 
edgeGain = 
            eg /. s; 
loopGain = 
            lg /. s; 
For[j = 1, 
            j ≤ n, \(j++\), 
z = 
              ev[\([j]\)]; 
rx = 
              Rn[\([j]\)]; 
lx = 
              Ln[\([j]\)]; 
fx = 
              lx . Bn . Un; 
gx = 
              Un . Cn . rx; 
linkEla = 
              linkElaMa /. s; 
loopEla = 
              cmrinv . linkEla; 
loopInf = 
              loopEla*z; 
For[k = 1, 
              k ≤ 
                Length[loopEla], \(k++\), 
AppendTo[
                outLoopEla, {
t, 
j, 
Re[z], 
Im[
                    z], 
k, 
loopsetText[\([k, 1]\)], 
loopnodes[\([k]\)], 
loopGain[\([k]\)], 
Re[
                    loopEla[\([k]\)]], 
Im[
                    loopEla[\([k]\)]], 
Abs[
                    loopEla[\([k]\)]], 
Re[
                    loopInf[\([k]\)]], 
Im[
                    loopInf[\([k]\)]], 
Abs[
                    loopInf[\([k]\)]]
}]
]
]
]; 
Export["outLinkEla.txt", 
          outLinkEla, "TSV"]; 
Export["outLoopEla.txt", outLoopEla, "TSV"];
];\)\)