LSStats: Difference between revisions

Latest revision as of 23:23, 8 November 2006

Summary

Screenshot

Download Model

Download the .scn, .sel and .lse files by clicking on the following link: File:LSStats.zip

Model Code Exploration

In the following sections we will examine all of the model files for this model. Note that instead of downloading the zip file above, you could just copy the text in the boxes below into a text editor and save it with the appropriate name (Section title). Opening the resulting .scn file in the SELES simulator would run this model.

stats.scn

SELES Scenario
$gisDir$ = ..\..\CaseStudy\gisData\cell
// load first layers
StudyArea = $baseGisDir$\StudyArea
// Open a patch layer
PatchLayer = $gisDir$\Old_bF
// Load model
Model Dimensions: PatchLayer
stats.sel
cwd oStats1
// Set up parameters
CellWidth = 100
HaPerCell = (CellWidth^2)/10000
NNType = 2 // For MPG
LandscapeId = 1 // This sets the variable LandscapeId to the integer represented by $x$
NumPTypes = 1
Reset Output $first$ // Tell SELES to not reset output files (i.e. to append to them) except on first iteration
// Run simulation
SimStart 10000 1

stats.sel

Seles Model
Landscape Events:
  identifyPatches.lse DEBUG
  stats.lse DEBUG
  centroid.lse DEBUG
  nn.lse DEBUG
Spatial Constants:
  PatchLayer
  StudyArea
Spatial Variables:
  EdgeInterior[2] <= 0
  PatchId[100000] <= 0
  CoreAreaId[100000] <= 0
  ComponentId[100000] <= 0
  NNLoc[-1] <= 0
  PatchLinks[3] <= 0
Global Constants:
 MaxPTypes = 4
 // Indices into patch stat vector
 NumPatchStats = 9
 rId = 0
 rType = 1
 rArea = 2
 rCoreArea = 3
 rPerim = 4
 rCentroid = 5
 rMaxCCE = 6
 rNNDist = 7
 rEffectiveId = 8
 // Indices into NN stat vector
 NumNNStats = 9 //8
 rEdgeType = 0
 rPatchType = 1
 rDist = 2
 rWeight = 3
 rStartLoc = 4
 rEndLoc = 5
 rNode1Id = 6
 rNode2Id = 7
 // Nearest neighbour algorithm types
 rNN = 0
 rMST = 1
 rMPG = 2
 MaxLoc = NumCells
 MaxDist = (NumRows^2 + NumCols^2)^(1/2) + 1
Global Variables:
 // PARAMETERS
 /* Cell size in ha and width in m */
 CellWidth = CELL WIDTH(PatchLayer)
 HaPerCell = (CellWidth * CellWidth) / 10000
 // Number of types to process (MUST BE LESS THAN MaxPTypes)
 /* For Planning Area, NumPTypes = 100 */
 /* For BEC, NumPTypes = 22 */
 /* For AgeClass, NumPTypes = 10 */
 /* For AgeClass2, NumPTypes = 4 */
//  NumPTypes = 4
 MinPType = 1
 NumPTypes = 1
 LandscapeId = 1
 Replicate = 1
 // Controls on nearest neighbour type
 NNType = rMPG
//  NNType = rMST
//  NNType = rNN
// Tracking variables
 A = 0
 NP = 0
 Area[MaxPTypes+1] = 0 OFF
 NumPatches[MaxPTypes+1] = 0 OFF
 MeanPS[MaxPTypes+1] = 0 OFF
 MaxPS[MaxPTypes+1] = 0 OFF
 MinPS[MaxPTypes+1] = 0 OFF
 PSSD[MaxPTypes+1] = 0 OFF
 TotalEdge[MaxPTypes+1] = 0 OFF
 TotalEdgeArea[MaxPTypes+1] = 0 OFF
 MSI[MaxPTypes+1] = 0 OFF
 AWMSI[MaxPTypes+1] = 0 OFF
 aan[MaxPTypes+1] = 0 OFF
 paRatio[MaxPTypes+1] = 0 OFF
 TCA[MaxPTypes+1] = 0 OFF
 NumCA[MaxPTypes+1] = 0 OFF
 MCA1[MaxPTypes+1] = 0 OFF
 CASD1[MaxPTypes+1] = 0 OFF
 MCA2[MaxPTypes+1] = 0 OFF
 CASD2[MaxPTypes+1] = 0 OFF
 MCAI[MaxPTypes+1] = 0 OFF
 IJI[MaxPTypes+1] = 0 OFF
 Eik[MaxPTypes+1, MaxPTypes+1] = 0 OFF
 NumNeighbs[MaxPTypes+1] = 0 OFF
 NumSameNeighbs[MaxPTypes+1] = 0 OFF
 // This should be in a separate section
 patchList[MaxPTypes+1] = 0 OFF
 coreAreaList[MaxPTypes+1] = 0 OFF
 //  patchList = 0 OFF
 pId = 0
 coreId = 0

identifyPatches.lse

/* Combination of stats from fragstats and apack */
/* Note: A particular problem should decide on stats required */
/* and adapt a custom stats event, which may include some of the */
/* stats herein, and may include new or modified stats. */
/* Apack stats are denoted using lower case identifiers. */
LSEVENT: STATS_LSE
DEFINITIONS
  LAYER: PatchLayer, StudyArea, PatchId, CoreAreaId
  LAYER: EdgeInterior
    FEATURE 1: Edge
    FEATURE 2: Interior   
  GLOBAL CONSTANT: MaxPTypes
  CONSTANT: pi = 3.141593
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId
     // Area indices
  GLOBAL VARIABLE: Area[]
  // Patch indices
  GLOBAL VARIABLE: NumPatches[]
  // Edge Indices
  GLOBAL VARIABLE: TotalEdge[], TotalEdgeArea[]
  // Core Area Indices
  GLOBAL VARIABLE: TCA[], NumCA[]
  // Adjacency matrix
  GLOBAL VARIABLE: Eik[]
  // pContagion
  GLOBAL VARIABLE: NumNeighbs[], NumSameNeighbs[]
  // Patch list
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rCoreArea, rPerim, rCentroid, rMaxCCE, rNNDist, MaxDist
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[], coreAreaList[]
  //   GLOBAL LIST{NumPatchStats} VARIABLE: patchList, coreAreaList
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  GLOBAL VARIABLE: pId, coreId
  EVENT VARIABLE: isPatchEvent, currType
  CLUSTER VARIABLE: numActiveCells, currPatchSize, currPatchCoreSize, currPatchPerim
  CLUSTER VARIABLE: rowTotal, colTotal
ENDDEF
INITIALSTATE
  EdgeInterior = 0
  pId = 0
  coreId = 0
  /* Three types of instances: the first NumPType instances identify and computes first order stats for patches */
  /*  (a separate instance is used per type so that the list is ordered by type) */
  /* The second identifies  and computes first order stats for core areas */
  /* The third computes second order stats and then summarizes and outputs */
  INITIALSTATE = 1 + NumPTypes
  isPatchEvent = EventId <= NumPTypes
  currType = EventId + MinPType - 1   // This is ok for patches, but don't use for core areas
ENDIS
RETURNTIME = IF Time EQ 0 THEN EventId + 365.25 ELSE 0
PROBINIT
  PROBINIT = IF (isPatchEvent) THEN (PatchLayer EQ currType) AND (StudyArea > 0)
             ELSE (EdgeInterior EQ Interior)
  numActiveCells = 0
  currPatchSize = 0
  currPatchCoreSize = 0
  currPatchPerim = 0
  rowTotal = 0
  colTotal = 0
  // For the patches ...
  IF (isPatchEvent)
     Area[PatchLayer] = Area[PatchLayer] + 1
     IF (PatchId EQ 0) // not visited...
        pId = pId + 1
        NumPatches[PatchLayer] = NumPatches[PatchLayer] + 1
     ENDFN
  // For the core areas ...
  ELSE
     coreId = coreId + (CoreAreaId EQ 0)
     NumCA[PatchLayer] = NumCA[PatchLayer] + (CoreAreaId EQ 0)
  ENDFN
ENDPI
TRANSITIONS
  TRANSITIONS = IF isPatchEvent THEN (PatchId EQ 0) ELSE (CoreAreaId EQ 0)
  // For the patches ...
  IF isPatchEvent
     numSimilarNeighbours = 0
     numDifferentNeighbours = 0
     OVER REGION CENTRED(1, 1)
        DECISION (0 <= PatchLayer < (MinPType + NumPTypes)) AND (StudyArea > 0)
        similarNeighb = (PatchLayer EQ currType)
        numSimilarNeighbours = numSimilarNeighbours + similarNeighb
        numDifferentNeighbours = numDifferentNeighbours + (similarNeighb EQ FALSE)
        NumNeighbs[PatchLayer] = NumNeighbs[PatchLayer] + 1
        NumSameNeighbs[PatchLayer] = NumSameNeighbs[PatchLayer] + similarNeighb 
        Eik[currType, PatchLayer] = Eik[currType, PatchLayer] + (similarNeighb EQ FALSE)
     ENDFN
     currPatchPerim = currPatchPerim + 4 - numSimilarNeighbours
     TotalEdge[PatchLayer] = TotalEdge[PatchLayer] + numDifferentNeighbours
     // Add in the diagonal neighbours
     OVER REGION CENTRED(1.1, 1.5)
        DECISION (0 <= PatchLayer < (MinPType + NumPTypes))
        similarNeighb = (PatchLayer EQ currType)
        numSimilarNeighbours = numSimilarNeighbours + similarNeighb 
        numDifferentNeighbours = numDifferentNeighbours + (similarNeighb EQ FALSE)
     ENDFN
     EdgeInterior = IF (numDifferentNeighbours > 0) THEN Edge ELSE Interior
     TCA[PatchLayer] = TCA[PatchLayer] + (EdgeInterior EQ Interior)
     TotalEdgeArea[PatchLayer] = TotalEdgeArea[PatchLayer] + (EdgeInterior EQ Edge)
     PatchId = pId
  ELSE
     CoreAreaId = coreId
  ENDFN
  currPatchSize = currPatchSize + 1
  currPatchCoreSize = currPatchCoreSize + (EdgeInterior EQ Interior)
  numActiveCells = numActiveCells + 1
  rowTotal = rowTotal + ROW(Location)
  colTotal = colTotal + COL(Location)
ENDTR
SPREADLOCATION
  REGION CENTRED(1, 1.5)
     DECISION IF isPatchEvent THEN (PatchId EQ 0) AND (PatchLayer EQ currType) AND (StudyArea > 0)
              ELSE (CoreAreaId EQ 0) AND (EdgeInterior EQ Interior) //AND (PatchLayer EQ currType)
ENDSL
SPREADTIMESTEP
  SPREADTIMESTEP = -1
  // If numActiveCells becomes 0, then this patch is done
  numActiveCells = numActiveCells - 1
  IF (numActiveCells EQ 0)
     // Add the patch to the list
     patchVar[rArea] = currPatchSize
     patchVar[rCoreArea] = currPatchCoreSize
     patchVar[rPerim] = currPatchPerim
     patchVar[rCentroid] = LOCATION(FLOOR(rowTotal/currPatchSize), FLOOR(colTotal/currPatchSize))
     patchVar[rMaxCCE] = 0
     patchVar[rNNDist] = -1
     IF isPatchEvent
        patchVar[rType] = currType
        patchVar[rId] = pId
        INSERT TAIL(patchList[currType], patchVar)
     ELSE
        patchVar[rType] = PatchLayer
        patchVar[rId] = coreId
        INSERT TAIL(coreAreaList[PatchLayer], patchVar)
     ENDFN
ENDFN
ENDST

stats.lse

/* Combination of stats from fragstats and apack */
/* Note: A particular problem should decide on stats required */
/* and adapt a custom stats event, which may include some of the */
/* stats herein, and may include new or modified stats. */
/* Apack stats are denoted using lower case identifiers. */
LSEVENT: STATS_LSE
DEFINITIONS
  LAYER: PatchLayer, StudyArea
  GLOBAL CONSTANT: MaxPTypes
  CONSTANT: pi = 3.141593
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId, Replicate
  /* Area indices */
  GLOBAL VARIABLE: A, Area[]
  /* Patch indices */
  GLOBAL VARIABLE: NP, MPS, tPSSD
  GLOBAL VARIABLE: NumPatches[], MeanPS[], MaxPS[], MinPS[], PSSD[]
  /* Edge Indices */
  GLOBAL VARIABLE: TotalEdge[], TotalEdgeArea[], MSI[], AWMSI[]
  /* related to inverse of MSI */
  GLOBAL VARIABLE: aan[], paRatio[]
  /* Core Area Indices */
  /* Total core area */
  GLOBAL VARIABLE: TCA[], NumCA[], MCA1[], CASD1[], MCA2[], CASD2[]
  GLOBAL VARIABLE: tMCA1, tCASD1, tMCA2, tCASD2, MCAI[]
  /* Interspersion and Juxtaposition */
  GLOBAL VARIABLE: IJI[]
  /* Adjacency matrix */
  GLOBAL VARIABLE: Eik[]
  /* pContagion */
  GLOBAL VARIABLE: NumNeighbs[], NumSameNeighbs[]
  // Patch list
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rCoreArea, rPerim, rCentroid, rMaxCCE, rNNDist, MaxDist
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[], coreAreaList[]
  //   GLOBAL LIST{NumPatchStats} VARIABLE: patchList, coreAreaList
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  MIN OUTPUT VARIABLE: LSStatsFile = "LSStats.txt"
  MIN OUTPUT VARIABLE: ClassStatsFile = "ClassStats.txt"
  //   MIN OUTPUT VARIABLE: PatchStatsFile = "PatchStats.txt"
ENDDEF
RETURNTIME
  RETURNTIME = IF Time EQ 0 THEN EventId + 1.5*365.25 ELSE 0
     // First compute second order stats for patches and core areas
     NP = 0
     A = 0
     MPS = 0
     tMCA1 = 0
     NCA = 0
     tMCA2 = 0
     OVER REGION WHOLE MAP
     // DECISION (0 <= PatchLayer < (MinPType + NumPTypes))
        DECISION StudyArea > 0
        A = A + 1
     ENDFN
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        NP = NP + NumPatches[pType]
     // A = A + Area[pType]
        MeanPS[pType] = IF (NumPatches[pType] > 0) THEN (Area[pType] / NumPatches[pType]) ELSE 0
        MPS = MPS + Area[pType]
        MCA1[pType] = IF (NumPatches[pType] > 0) THEN (TCA[pType] / NumPatches[pType]) ELSE 0
        tMCA1 = tMCA1 + TCA[pType]
        NCA = NCA + NumCA[pType]
        MCA2[pType] = IF (NumCA[pType] > 0) THEN (TCA[pType] / NumCA[pType]) ELSE 0
        tMCA2 = tMCA2 + TCA[pType]
     ENDFN
     MPS = MPS / NP
     tMCA1 = tMCA1 / NP
     tMCA2 = tMCA2 / NCA
     // Iterate over list of patches to compute stats
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        i = Index + MinPType
        pos = HEAD(patchList[i])
        IF pos
           MaxPS[i] = patchVar[rArea]
           MinPS[i] = patchVar[rArea]
        ENDFN
        WHILE (pos > 0)
           patchVar [=] GET(patchList[i], pos)
           pos = NEXT(patchList[i], pos)
     //  i = patchVar[rType]
           patchSize = patchVar[rArea]
           patchCoreSize = patchVar[rCoreArea]
           patchPerim = patchVar[rPerim]
           MaxPS[i] = MAX(MaxPS[i], patchSize)
           MinPS[i] = MIN(MinPS[i], patchSize)
     // Sum square of different between patch size and mean patch size
           PSSD[i] = PSSD[i] + (patchSize - MeanPS[i])^2
           tPSSD = tPSSD + (patchSize - MPS)^2
           MSI[i] = MSI[i] + (0.25 * patchPerim / (patchSize ^ (1/2)))
           AWMSI[i] = AWMSI[i] + (0.25 * patchPerim * (patchSize ^ (1/2)))
           aan[i] = aan[i] + 16 * patchSize  / (patchPerim^2)
           paRatio[i] = paRatio[i] + patchPerim / patchSize
           CASD1[i] = CASD1[i] + (patchCoreSize - MCA1[i])^2
           tCASD1 = tCASD1 + (patchCoreSize - tMCA1)^2
           MCAI[i] = IF (patchSize > 0) THEN (MCAI[i] + patchCoreSize / patchSize) ELSE MCAI[i]
     /*    OUTPUT RECORD(PatchStatsFile)
              LandscapeId: LandscapeId
              Replicate: Replicate
              id: patchVar[rId]
              type: patchVar[rType]
              size: patchVar[rArea]
              coreSize: patchVar[rCoreArea]
           ENDFN
     */
        ENDFN
        // Iterate over list of core areas to compute stats
        /* Sum square of different between patch size and mean patch size */
        pos = HEAD(coreAreaList[i])
        WHILE (pos)
           patchVar [=] GET(coreAreaList[i], pos)
           pos = NEXT(coreAreaList[i], pos)
        // i = patchVar[rType]
           patchSize = patchVar[rArea]
           patchCoreSize = patchVar[rCoreArea]
           patchPerim = patchVar[rPerim]
           CASD2[i] = CASD2[i] + (patchCoreSize - MCA2[i])^2
           tCASD2 = tCASD2 + (patchCoreSize - tMCA2)^2
        ENDFN
        ENDFN
     // Finally summarize over patch types and output results
     tTE = 0
     tTEA = 0
     MaxPatch = 0
     MinPatch = NUMCELLS
     tMSI = 0 
     tAWMSI = 0 
     tMCAI = 0
     taan = 0
     taan2 = 0
     tap2 = 0
     tTCA = 0
     NCA = 0
     tIJI = 0
     tContag = 0
     tContag2 = 0
     measuredDiversity = 0
     measuredDiversity2 = 0
     do = 0
     SHDI = 0
     SIDI = 1
     MSIDI = 0
     PR = 0
     asm = 0
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        tTE = tTE + TotalEdge[pType]
        tTEA = tTEA + TotalEdgeArea[pType]
        tap2 = IF (NumPatches[pType] > 0) THEN tap2 + TotalEdge[pType] / NumPatches[pType] ELSE tap2
        MaxPatch = MAX(MaxPatch, MaxPS[pType])
        MinPatch = MIN(MinPatch, MinPS[pType])
        PSSD[pType] = IF (NumPatches[pType] > 0) THEN (PSSD[pType] / NumPatches[pType])^(1/2) ELSE 0
        tMSI = tMSI + MSI[pType]
        MSI[pType] = IF (NumPatches[pType] > 0) THEN MSI[pType] / NumPatches[pType] ELSE 0
        tAWMSI = tAWMSI + AWMSI[pType]
        AWMSI[pType] = IF (Area[pType] > 0) THEN AWMSI[pType] / Area[pType] ELSE 0
        taan  = taan  + aan[pType]
        aan[pType] = IF (NumPatches[pType] > 0) THEN aan[pType] / NumPatches[pType] ELSE 0
        taan2  = taan2  + aan[pType]
        paRatio[pType] = IF (NumPatches[pType] > 0) THEN paRatio[pType] / NumPatches[pType] ELSE 0
        tTCA = tTCA + TCA[pType]
        NCA = NCA + NumCA[pType]
        CASD1[pType] = IF (NumPatches[pType] > 0) THEN (CASD1[pType] / NumPatches[pType])^(1/2) ELSE 0
        CASD2[pType] = IF (NumCA[pType] > 0) THEN (CASD2[pType] / NumCA[pType])^(1/2) ELSE 0 
        tMCAI = tMCAI + MCAI[pType]
        PR = PR + (Area[pType] NEQ 0)
     ENDFN
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        OVER INDEX SEQUENCE(pType + 1, NumPTypes - 1)
           pType2 = Index
           x = IF (tTE > 0) THEN Eik[pType, pType2] / tTE ELSE 0
           tIJI = tIJI + x * LOG(x)
        ENDFN
        Pi = Area[pType] / A
        OVER INDEX SEQUENCE(0, NumPTypes-1)
           pType2 = Index + MinPType
           AMik = IF (TotalEdge[pType] > 0) THEN Eik[pType, pType2] / TotalEdge[pType] ELSE 0
           IJI[pType] = IJI[pType] + AMik * LOG(AMik)
           x = IF ((TotalEdge[pType] > 0) AND (Pi > 0)) THEN Pi * Eik[pType, pType2] / TotalEdge[pType] ELSE 0
           tContag = tContag + x * LOG(x)
           measuredDiversity  = measuredDiversity - AMik * LOG(AMik)
           measuredDiversity2  = IF(pType NEQ pType2) THEN  measuredDiversity2 - AMik * LOG(AMik) ELSE measuredDiversity2
           asm = asm + AMik^2
        ENDFN
        IJI[pType] = IF (PR > 2) THEN (-1 * IJI[pType]/LOG(PR-1)) ELSE 0
        SHDI = IF (Area[pType] EQ 0) THEN SHDI ELSE SHDI - (Area[pType]/A)*LOG(Area[pType]/A)
        SIDI = SIDI - (Area[pType]/A)^2
        MSIDI = IF (Area[pType] EQ 0) THEN MSIDI ELSE MSIDI - LOG((Area[pType]/A)^2)
     ENDFN
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        Pi = Area[pType] / A
        /* Assume that Pi < 0.5 */
        Contag2 = (NumSameNeighbs[pType] / NumNeighbs[pType]) - Pi
        Contag2 = IF (Pi <= 0) OR (Pi >= 1) THEN 0
                  ELSE IF (Contag2 >= 0) THEN Contag2 / (1 - Pi) ELSE Contag2 / Pi
        do = do - Pi * LOG(Pi) 
        OUTPUT RECORD(ClassStatsFile)
           DECISION Area[pType] > 0
           LandscapeId: LandscapeId
           Replicate: Replicate
           pType: pType
           A: Area[pType] * HaPerCell
           PCTLAND: 100 * Pi
           LPI: 100 * MaxPS[pType] / A
           LargestPatch: MaxPS[pType] * HaPerCell
           SmallestPatch: MinPS[pType] * HaPerCell
           NP: NumPatches[pType]
           PD: 100 * NumPatches[pType] / (A * HaPerCell)
           MPS: MeanPS[pType] * HaPerCell
           PSSD: PSSD[pType] * HaPerCell
           PSCV: 100 * PSSD[pType]/MeanPS[pType]
           TE: CellWidth * TotalEdge[pType]
           ED: CellWidth * TotalEdge[pType] / A
           ap: CellWidth * TotalEdge[pType] / NumPatches[pType]
           TEA: TotalEdgeArea[pType] * HaPerCell
           EDA: 100 * TotalEdgeArea[pType] / A
           OPOE: TotalEdgeArea[pType] / (TotalEdge[pType] * NumPatches[pType])
           LSI: 0.25 * TotalEdge[pType] / ((Area[pType])^(1/2))
           MSI: MSI[pType]
           AWMSI: AWMSI[pType]
           aan: aan[pType]
           CPCTLAND: 100 * TCA[pType] / A
           TCA: TCA[pType] * HaPerCell
           NCA: NumCA[pType]
           CAD: 100 * NumCA[pType] / (A * HaPerCell)
           MCA1: MCA1[pType] * HaPerCell
           CASD1: CASD1[pType] * HaPerCell
           CACV1: 100 * CASD1[pType]/MCA1[pType]
           MCA2: MCA2[pType] * HaPerCell
           CASD2: CASD2[pType] * HaPerCell
           CACV2: 100 * CASD2[pType]/MCA2[pType]
           TCAI: 100 * TCA[pType] / Area[pType]
           MCAI: 100 * MCAI[pType] / NumPatches[pType]
           IJI: 100 * IJI[pType]
           /* Zero -> random adjacency (in [-1, 1] */
           /* Positive -> more clumped than random */
           /* Negative -> less clumped than random */
           CONTAG2: Contag2
           NumNeighbs: NumNeighbs[pType]
           NumSameNeighbs: NumSameNeighbs[pType]
        ENDFN
        tContag2 = tContag2 + Contag2
     ENDFN
     tPSSD = (tPSSD / NP)^(1/2)
     tCASD1 = (tCASD1 / NP)^(1/2)
     tCASD2 = (tCASD2 / NCA)^(1/2)
     tIJI = IF (PR > 2) THEN (-1 * tIJI/LOG((1/2) * PR * (PR-1))) ELSE 0
     tContag = 1 + tContag / (2 * LOG(PR))
     OUTPUT RECORD(LSStatsFile)
        LandscapeId: LandscapeId
        Replicate: Replicate
        A: A * HaPerCell
        LPI: 100 * MaxPatch / A
        LargestPatch: MaxPatch  * HaPerCell
        SmallestPatch: MinPatch  * HaPerCell
        NP: NP
        PD: 100 * NP  / (A * HaPerCell)
        MPS: MPS * HaPerCell
        PSSD: tPSSD * HaPerCell
        PSCV: 100 * tPSSD/MPS
        TE: CellWidth * tTE
        ED: CellWidth * tTE / A
        ap: CellWidth * tTE / NP
        ap2: CellWidth * tap2
        TEA: tTEA * HaPerCell
        EDA: 100 * tTEA / A
        OPOE: tTEA / (tTE * NP)
        LSI: 0.25 * tTE / (A^(1/2))
        MSI: tMSI / NP
        AWMSI: tAWMSI / A
        aan: taan / NP
        aan2: taan2
        TCA: tTCA * HaPerCell
        NCA: NCA
        CAD: 100 * NCA / (A * HaPerCell)
        MCA1: tMCA1 * HaPerCell
        CASD1: tCASD1 * HaPerCell
        CACV1: 100 * tCASD1/tMCA1
        MCA2: tMCA2 * HaPerCell
        CASD2: tCASD2 * HaPerCell
        CACV2: 100 * tCASD2/tMCA2
        TCAI: 100 * tTCA / A
        MCAI: 100 * tMCAI / NP
        tIJI: 100 * tIJI
        CONTAG1: 100 * tContag
        CONTAG2: tContag2 / PR
        co: 2 * LOG(PR) - measuredDiversity
        co2: (LOG(PR^2 + PR) - LOG(2)) - measuredDiversity
        col1: 1 - measuredDiversity / (2*LOG(PR))
        cor: 1 - measuredDiversity / (LOG(PR^2 + PR) - LOG(2))
        asm: asm
        do: LOG(PR) - do
        dor: 1 - LOG(PR) / do
        ede: 1 - measuredDiversity2 / (2*LOG(PR))
        ede2: 1 - measuredDiversity2 / (LOG(PR^2 + PR) - LOG(2))
        SHDI: SHDI
        SIDI: SIDI
        MSIDI: MSIDI
        PR: PR
        PRD: 100 * PR / (A * HaPerCell)
        RPR: PR / (NumPTypes - 1)
        SHEI: IF (PR > 1) THEN SHDI / LOG(PR) ELSE 0
        SIEI: IF (PR > 1) THEN SIDI / (1 - 1/PR) ELSE 0
        MSIEI: IF (PR > 1) THEN MSIDI / LOG(PR) ELSE 0
     ENDFN
ENDRT
// Don't continue
NUMCLUSTERS = 0

centroid.lse

LSEVENT: CENTROID_LSE
DEFINITIONS
  GLOBAL CONSTANT: MaxPTypes
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId, Replicate
  /* Assume these stats are defined */
  GLOBAL VARIABLE: A, Area[], NP, NumPatches[]
  /* Mean Connectivity between patches */
  GLOBAL VARIABLE: MeanCCE[MaxPTypes+1], MaxCCE[MaxPTypes+1], MeanCD[MaxPTypes+1]
  GLOBAL VARIABLE: tCCE, tMaxCCE, tCD
  // Patch list: type and area previously defined
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rCentroid, rMaxCCE
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[]
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  EVENT VARIABLE: id
  CLUSTER VARIABLE: rowTotal, colTotal, numActiveCells
//   MIN OUTPUT VARIABLE: LSStatsFile = "LSStatsC.txt"
  MIN OUTPUT VARIABLE: ClassStatsFile = "ClassStatsC.txt"
//   MIN OUTPUT VARIABLE: EdgeStatsFile = "EdgeStatsC.txt"
ENDDEF
INITIALSTATE
  tCCE = 0
  tMaxCCE = 0
  tCD = 0
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     pType = Index + MinPType
     MeanCCE[pType] = 0
     MaxCCE[pType] = 0
     MeanCD[pType] = 0
  ENDFN
  /* All the stats area available in the patch list, we only need to */
  /* summarize and output results */
  INITIALSTATE = 1
ENDIS
RETURNTIME
//   RETURNTIME = IF Time EQ 0 THEN (4 * 365.25) ELSE 0
  RETURNTIME = IF Time EQ 0 THEN (4 * 365.25) ELSE 0
  /* Finalize stats */
  // Iterate over each pair of patches of the same type
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     i = Index + MinPType
     pos1 = HEAD(patchList[i])
     WHILE (pos1)
        patchVar [=] GET(patchList[i], pos1)
        currType = patchVar[rType]
        id1 = patchVar[rId]
        patchArea1 = patchVar[rArea]
        centroid1 = patchVar[rCentroid]
        maxCCE1 = patchVar[rMaxCCE]
        nextPos = NEXT(patchList[i], pos1)
//         pos2 = HEAD(patchList[i])
        pos2 = nextPos
        WHILE (pos2)
           IF (currType EQ GET(patchList[i], pos2, rType))
              patchVar [=] GET(patchList[i], pos2)
              patchArea2 = patchVar[rArea]
              centroid2 = patchVar[rCentroid]
              maxCCE2 = patchVar[rMaxCCE]
             /* Determine the max patch CCE */
              d = DISTANCE(centroid1, centroid2)
              CCE = (patchArea1 * patchArea2) / d^2
              MeanCCE[currType] = MeanCCE[currType] + CCE
              maxCCE1 = MAX(CCE, maxCCE1)
              maxCCE2 = MAX(CCE, maxCCE2)
              MeanCD[currType] = MeanCD[currType] + d
              // Update maxCCE in patches
              SET(patchList[i], pos2, rMaxCCE, maxCCE2)
/*
              OUTPUT RECORD(EdgeStatsFile)
                 LandscapeId: LandscapeId
                 Replicate: Replicate
                 type: currType
                 patchId1: id1
                 patchId2: patchVar[rId]
                 d: d
              ENDFN
*/
           ENDFN
           pos2 = NEXT(patchList[i], pos2)
        ENDFN
       /* At this point, the patch CCE for p1 will be known */
        SET(patchList[i], pos1, rMaxCCE, maxCCE1)
        MaxCCE[currType] = MaxCCE[currType] + maxCCE1
        pos1 = nextPos
     ENDFN
  ENDFN
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     currType = Index + MinPType
     tCCE = tCCE + MeanCCE[currType]
     tMaxCCE = tMaxCCE + MaxCCE[currType]
     MeanCCE[currType] = MeanCCE[currType] / (NumPatches[currType] * (NumPatches[currType] - 1) * 0.5)
     MaxCCE[currType] = MaxCCE[currType] / (NumPatches[currType] - 1)
     tCD = tCD + MeanCD[currType]
     MeanCD[currType] = MeanCD[currType] / (NumPatches[currType] * (NumPatches[currType] - 1) * 0.5)
     OUTPUT RECORD(ClassStatsFile)
        DECISION Area[currType] > 0
        LandscapeId: LandscapeId
        Replicate: Replicate
        pType: currType
        CCE: MeanCCE[currType]
        MaxCCE: MaxCCE[currType]
        CD: MeanCD[currType]
     ENDFN
  ENDFN
  tCCE = IF (NP <= 1) THEN 0 ELSE tCCE / (NP * (NP - 1) * 0.5)
  tMaxCCE = tMaxCCE / NP
  tCD = IF (NP <= 1) THEN 0 ELSE tCD / (NP * (NP - 1) * 0.5)
/*
  OUTPUT RECORD(LSStatsFile)
     LandscapeId: LandscapeId
     Replicate: Replicate
     tCCE: tCCE
     tMaxCCE: tMaxCCE
     tCD: tCD
  ENDFN
*/
ENDRT
EVENTOPENINGS = 0

nn.lse

LSEVENT: NearestNeighbour_LSE
DEFINITIONS
  // Nearest neighbour algorithm types
  GLOBAL CONSTANT: rNN, rMST, rMPG
  GLOBAL VARIABLE: NNType
  /* Assume these layers are all defined */
  LAYER: PatchLayer, PatchId, StudyArea
  LAYER: EdgeInterior
    FEATURE 1: Edge
    FEATURE 2: Interior
  // Scratch id layer for connectve components
  LAYER: ComponentId, PatchLinks 
  GLOBAL CONSTANT: MaxPTypes, MaxDist
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId, Replicate
  LAYER: NNLoc
  GLOBAL CONSTANT: MaxLoc
  /* Assume these stats are defined */
  GLOBAL VARIABLE: A, Area[], NP, NumPatches[]
  // Patch list: previously defined (except rNNDist, which is at MaxDist)
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rNNDist, rEffectiveId
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[]
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  /* Nearest Neighbour */
  GLOBAL VARIABLE: MeanNN[MaxPTypes+1], MinNN[MaxPTypes+1], MaxNN[MaxPTypes+1], NNSD[MaxPTypes+1]
  GLOBAL VARIABLE: MeanNN2[MaxPTypes+1]
  GLOBAL VARIABLE: tMeanNN, tMinNN, tMaxNN
  /* MST */
  GLOBAL VARIABLE: MST[MaxPTypes+1], meanMST[MaxPTypes+1]
  GLOBAL VARIABLE: startId[MaxPTypes+1]
  /* Straigh line least cost */
  GLOBAL VARIABLE: MPG[MaxPTypes+1], meanMPG[MaxPTypes+1], nMPG[MaxPTypes+1]
  GLOBAL VARIABLE: AW_MPG[MaxPTypes+1], meanAW_MPG[MaxPTypes+1], nAW_MPG[MaxPTypes+1]
  /* NOTE: This variable must be resized to NP (which is unkown before running stats.lse) */
  //   GLOBAL VARIABLE: effectivePatchId[2], PosList[2]
  GLOBAL VARIABLE: effectivePatchId[2], PosList[200]
  // NN graph 
  GLOBAL CONSTANT: NumNNStats, rEdgeType, rPatchType, rDist, rWeight, rStartLoc, rEndLoc, rNode1Id, rNode2Id
  GLOBAL GRAPH{NumPatchStats, NumNNStats} VARIABLE: nnGraph[MaxPTypes+1]
  GLOBAL VARIABLE: edgeVar[NumNNStats], edgeVar2[NumNNStats]
  EVENT VARIABLE: currType, PatchesToConnect, ComponentsToConnect
  CLUSTER VARIABLE: currPatchId, patchPos
  CELL VARIABLE: currEffectivePid, newEffectivePid
  CELL VARIABLE: startLocation, numSteps
  //   MIN OUTPUT VARIABLE: LSStatsFile = "LSStatsNN.txt"
  MIN OUTPUT VARIABLE: ClassStatsFile = "ClassStatsNN.txt"
  //   MIN OUTPUT VARIABLE: PatchStatsFile = "PatchStatsNN.txt"
  //   MIN OUTPUT VARIABLE: EdgeStatsFile = "EdgeStatsNN.txt"
ENDDEF
INITIALSTATE
  /* The first NumPType instances compute stats for each type */
  /* The last summarizes and outputs */
  INITIALSTATE = NumPTypes + 1
  currType = EventId + MinPType - 1
ENDIS
RETURNTIME
  //   RETURNTIME = IF Time EQ 0 THEN (EventId + 5 * 365.25) ELSE 0
  RETURNTIME = IF Time EQ 0 THEN (EventId + 6 * 365.25) ELSE 0
  /* Finalize and output stats */
  IF currType >= NumPTypes + MinPType
     tMeanNN = 0
     tMinNN = MaxDist
     tMaxNN = 0
     N_MST = 0
     tMeanMST = 0
     tMST = 0
     tAWMST = 0
     tMPG = 0
     tMeanMPG = 0
    // Iterate over list of patches to compute stats
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        MeanNN[pType] = 0
        MinNN[pType] = MaxDist
        MaxNN[pType] = 0
        NNSD[pType] = 0
        MeanNN2[pType] = 0
        nNN = 0
        MST[pType] = 0
meanMST[pType] = 0
        MPG[pType] = 0
        meanMPG[pType] = 0
        nMPG[pType] = 0
        AW_MPG[pType] = 0
        meanAW_MPG[pType] = 0
        nAW_MPG[pType] = 0
        pos = FIRST(nnGraph[pType])
        WHILE (pos)
           patchVar [=] GET(nnGraph[pType], pos)
           pos = NEXT(nnGraph[pType], pos)
     //            pType = patchVar[rType]
           MeanNN[pType] = MeanNN[pType] + patchVar[rNNDist]
           IF (patchVar[rNNDist] > 0)
              MinNN[pType] = MIN(patchVar[rNNDist], MinNN[pType])
              MaxNN[pType] = MAX(patchVar[rNNDist], MaxNN[pType])
           ENDFN
     /*
           OUTPUT RECORD(PatchStatsFile)
              LandscapeId: LandscapeId
              Replicate: Replicate
              id: patchVar[rId]
              patchType: pType
              nnDist: patchVar[rNNDist] * CellWidth
           ENDFN
     */
        ENDFN
        // Iterate over list of edges to compute stats
        pos = FIRST LINK(nnGraph[pType])
        WHILE (pos)
           edgeVar [=] GET LINK(nnGraph[pType], pos)
           pos = NEXT LINK(nnGraph[pType], pos)
     //            pType = edgeVar[rPatchType]
IF (edgeVar[rEdgeType] EQ rNN)
 MeanNN2[pType] = MeanNN2[pType] + edgeVar[rDist]
 nNN = nNN + 1
ENDFN
           IF (edgeVar[rEdgeType] EQ rNN) OR (edgeVar[rEdgeType] EQ rMST)
              MST[pType] = MST[pType] + edgeVar[rDist]
           ENDFN
           MPG[pType] = MPG[pType] + edgeVar[rDist]
           nMPG[pType] = nMPG[pType] + 1
           edgeWeight = (CellWidth * edgeVar[rDist]) / (HaPerCell * edgeVar[rWeight])  // m/ha
           AW_MPG[pType] = AW_MPG[pType] + edgeWeight
           nAW_MPG[pType] = nAW_MPG[pType] + edgeWeight
     /*
           OUTPUT RECORD(EdgeStatsFile)
              LandscapeId: LandscapeId
              Replicate: Replicate
              patchType: pType
              edgeType: edgeVar[rEdgeType]
              id1: edgeVar[rNode1Id]
              id2: edgeVar[rNode2Id]
              dist: edgeVar[rDist] * CellWidth
              edgeWeight: edgeWeight
              node1Row: ROW(edgeVar[rStartLoc])
              node1Col: COL(edgeVar[rStartLoc])
              node2Row: ROW(edgeVar[rEndLoc])
              node3Col: COL(edgeVar[rEndLoc])
           ENDFN
     */
        ENDFN
        tMeanNN = tMeanNN + MeanNN[pType]
        IF MinNN[pType] > 0
           tMinNN = MIN(tMinNN, MinNN[pType])
        ENDFN
        tMaxNN = MAX(tMaxNN, MaxNN[pType])
        IF NumPatches[pType] > 0
           MeanNN[pType] = MeanNN[pType] / NumPatches[pType]
        ENDFN
        IF nNN > 0
           MeanNN2[pType] = MeanNN2[pType] / nNN
        ENDFN
        IF (NumPatches[pType] > 1)
           meanMST[pType] = MST[pType] / (NumPatches[pType] - 1)
           tMeanMST = tMeanMST + meanMST[pType]
           tMST = tMST + MST[pType]
           tAWMST = tAWMST + MST[pType] * Area[pType]
           N_MST = N_MST + 1
        ENDFN
        IF (nMPG[pType] > 1)
           meanMPG[pType] = MPG[pType] / nMPG[pType]
           tMPG = tMPG + MPG[pType]
        ENDFN
        IF (nAW_MPG[pType] > 1)
           meanAW_MPG[pType] = AW_MPG[pType] / nAW_MPG[pType]
           tAW_MPG = tAW_MPG + AW_MPG[pType]
        ENDFN
     ENDFN
     tMST = IF (N_MST > 0) THEN tMST / N_MST ELSE 0
     tAWMST = IF (A > 0) THEN tAWMST / A ELSE 0
     tMPG = IF (N_MST > 0) THEN tMPG / N_MST ELSE 0
     /* Compute stddev and second order stats */
     tNNSD = 0
     tMSTSD = 0
     tMPGSD = 0
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
   /*
        pos = FIRST LINK(nnGraph[pType])
        WHILE (pos)
           edgeVar [=] GET LINK(nnGraph[pType], pos)
           pos = NEXT LINK(nnGraph[pType], pos)
           IF (edgeVar[rEdgeType] EQ rNN)
              NNSD[pType] = NNSD[pType] + (edgeVar[rDist] - MeanNN[pType])^2
           ENDFN
        ENDFN
  */
        pos = FIRST(nnGraph[pType])
        WHILE (pos)
           patchVar [=] GET(nnGraph[pType], pos)
           pos = NEXT(nnGraph[pType], pos)
//            pType = patchVar[rType]
           NNSD[pType] = NNSD[pType] + (patchVar[rNNDist] - MeanNN[pType])^2
        ENDFN
        tNNSD = tNNSD + NNSD[pType]
        NNSD[pType] = IF (NumPatches[pType] EQ 0) THEN 0 ELSE (NNSD[pType] / NumPatches[pType])^(1/2)
        tMSTSD  = tMSTSD  + (MST[pType] - tMST)^2
        tMPGSD  = tMPGSD  + (MPG[pType] - tMPG)^2
        OUTPUT RECORD(ClassStatsFile)
           DECISION Area[pType] > 0
           LandscapeId: LandscapeId
           Replicate: Replicate
           pType: pType
           MNN: CellWidth * MeanNN[pType]
           MNN2: CellWidth * MeanNN2[pType]
           MinNN: CellWidth * MinNN[pType]
           MaxNN: CellWidth * MaxNN[pType]
           NNSD: CellWidth * NNSD[pType]
           NNCV: 100 * NNSD[pType] / MeanNN[pType]
           Dispersion: 2 * (NumPatches[pType] / Area[pType])^(1/2) * MeanNN[pType]
           meanMST: IF (NNType EQ rNN) THEN 0 ELSE CellWidth * meanMST[pType]
           tMST: IF (NNType EQ rNN) THEN 0 ELSE CellWidth * MST[pType]
           meanMPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * meanMPG[pType]
           tMPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * MPG[pType]
           nMPG: nMPG[pType]
           meanAW_MPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * meanAW_MPG[pType]
           tAW_MPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * AW_MPG[pType]
           nAW_MPG: nMPG[pType]
        ENDFN
     ENDFN
     tMeanNN = tMeanNN / NP
     tNNSD = (tNNSD  / NP)^(1/2)
     tMSTSD  = (tMSTSD  / N_MST)^(1/2)
     tMPGSD  = (tMPGSD  / N_MST)^(1/2)
/*
     OUTPUT RECORD(LSStatsFile)
        LandscapeId: LandscapeId
        Replicate: Replicate
        MNN: CellWidth * tMeanNN
        MinNN: CellWidth * tMinNN
        MaxNN: CellWidth * tMaxNN
        NNSD: CellWidth * tNNSD 
        NNCV: 100 * tNNSD / tMeanNN
        MST: IF (NNType EQ rNN) THEN 0 ELSE CellWidth * tMST
//         tMSTSD: IF (NNType EQ rNN) THEN 0 ELSE CellWidth * tMSTSD
//         MSTCV: IF (NNType EQ rNN) THEN 0 ELSE 100 * tMSTSD / tMST
        AWMST: IF (NNType EQ rNN) THEN 0 ELSE CellWidth * tAWMST
        MPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * tMPG
        AWMPG: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * tAW_MPG
//         tMPGSD: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE CellWidth * tMPGSD
//         MPGCV: IF (NNType EQ rNN) OR (NNType EQ rMST) THEN 0 ELSE 100 * tMPGSD / tMPG
     ENDFN
*/
  ENDFN
  IF ((currType < (NumPTypes + MinPType)) AND (NumPatches[currType] > 1))
///      PatchId = 0
     ComponentId = 0
     NNLoc = MaxLoc
     PatchesToConnect = NumPatches[currType]
     ComponentsToConnect = NumPatches[currType] - 1
  ENDFN
  /* NOTE: BE very careful to be sure that variables shared with */
  /* other events are valid (especially after initialization in INITIALSTATE!) */
  IF (currType EQ MinPType)
     /* Resize the effective patch id vector */
     RESIZE(effectivePatchId, NP+1)
     RESIZE(PosList, NP+1)
//RESIZE(effectivePatchId, 200)
//RESIZE(PosList, 200)
     tmpPid = 1
     OVER INDEX SEQUENCE(0, NumPTypes-1)
        pType = Index + MinPType
        OVER INDEX SEQUENCE(0, NumPatches[pType]-1)
           effectivePatchId[tmpPid] = tmpPid
           tmpPid = tmpPid + 1
        ENDFN
        // NOTE: the patches must be ordered by type for this to work
        startId[pType] = IF (pType EQ 0) THEN 0 ELSE startId[pType-1] + NumPatches[pType-1]
/*
        // Copy nodes to graph.  Alternative: define patchList as a graph
        pos = FIRST(patchList[pType])
        WHILE (pos)
           patchVar [=] GET(patchList[pType], pos)
           pos = NEXT(patchList[pType], pos)
           patchVar[rNNDist] = MaxDist
           INSERT(nnGraph[pType], patchVar)
           i = patchVar[rId]
           PosList[i] = FIRST(nnGraph[pType])
        ENDFN
*/
     ENDFN
  ENDFN
ENDRT
/* Don't continue if this is the output instance */
NUMCLUSTERS
//   IF (currType EQ NumPTypes) OR (NumPatches[currType] <= 1) THEN 0 ELSE -1
  NUMCLUSTERS = IF (currType >= (NumPTypes+MinPType)) OR (NumPatches[currType] <= 1) THEN 0 ELSE -1
ENDNC
PROBINIT
//   PROBINIT = (PatchLayer EQ currType) AND (EdgeInterior EQ Edge)
  PROBINIT = (PatchLayer EQ currType) AND (StudyArea > 0)
  currPatchId = PatchId
  startLocation = Location
  numSteps = 0
  pos = FIND(patchList[currType], patchVar, patchVar[rId] EQ currPatchId)
  pos2 = FIND(nnGraph[currType], patchVar, patchVar[rId] EQ currPatchId)
  patchVar [=] GET(patchList[currType], pos)
  i = patchVar[rId]
  IF (!pos2)
     patchVar[rNNDist] = MaxDist
     INSERT(nnGraph[currType], patchVar)
     PosList[i] = FIRST(nnGraph[currType])
  ENDFN
  patchPos = PosList[i]
///   patchPos = FIND(nnGraph[currType], patchVar, patchVar[rId] EQ currPatchId)
//   patchPos = PosList[currPatchId]
ENDPI
TRANSITIONS
  /* Check if we meet another spreading component */
  IF ((ComponentId NEQ 0) AND (effectivePatchId[ComponentId] NEQ effectivePatchId[currPatchId]))
     // Find the patch nodes
//      nPos2 = FIND(nnGraph[currType], patchVar, patchVar[rId] EQ ComponentId)
     nPos2 = PosList[ComponentId]
//      d = DISTANCE(NNLoc, startLocation)
     rowDiff = MAX(0, |ROW(NNLoc) - ROW(startLocation)| - 1)
     colDiff = MAX(0, |COL(NNLoc) - COL(startLocation)| - 1)
     d = (rowDiff^2 + colDiff^2)^0.5
          // Add the patch to the list or update the distance if it is already there
//      edgePos = FIND LINK(nnGraph[currType], edgeVar, ((edgeVar[rNode1Id] EQ currPatchId) AND (edgeVar[rNode2Id]  EQ ComponentId)) OR ((edgeVar[rNode1Id] EQ ComponentId) AND (edgeVar[rNode2Id] EQ currPatchId)))
     // Set NN at patch level
     currD1 = GET(nnGraph[currType], patchPos, rNNDist)
     IF d < currD1
        SET(nnGraph[currType], patchPos, rNNDist, d)
     ENDFN
     currD2 = GET(nnGraph[currType], nPos2, rNNDist)
     IF d < currD2
        SET(nnGraph[currType], nPos2, rNNDist, d)
     ENDFN
     PatchesToConnect = PatchesToConnect - (currD1 EQ MaxDist) - (currD2 EQ MaxDist)
     // Create a link
     edgeVar[rEdgeType] = IF (currD1 EQ MaxDist) OR (currD2 EQ MaxDist) THEN rNN ELSE rMST
     edgeVar[rPatchType] = PatchLayer
     edgeVar[rNode1Id] = currPatchId
     edgeVar[rNode2Id] = ComponentId
     edgeVar[rDist] = d
     Area1 = GET(nnGraph[currType], patchPos, rArea)
     Area2 = GET(nnGraph[currType], nPos2, rArea)
     edgeVar[rWeight] = (Area1 * Area2)^0.5
     edgeVar[rStartLoc] = startLocation
     edgeVar[rEndLoc] = NNLoc
     INSERT LINK(nnGraph[currType], patchPos, nPos2, edgeVar)
     IF (currType EQ MinPType)
        OVER REGION VECTOR(startLocation,NNLoc)
           PatchLinks = edgeVar[rEdgeType]+1
        ENDFN
     ENDFN
     ComponentsToConnect = ComponentsToConnect - 1
     /* Need to update all effective id's */
     currEffectivePid = effectivePatchId[currPatchId]
     newEffectivePid = effectivePatchId[ComponentId]
     OVER INDEX SEQUENCE(1, NP)
        // NOTE: the patches must be ordered by type for this to work
//      OVER INDEX SEQUENCE(startId[currType], startId[currType] + NumPatches[currType] - 1)
        effectivePatchId[Index] = IF (effectivePatchId[Index] EQ currEffectivePid)
                                  THEN newEffectivePid ELSE effectivePatchId[Index]
     ENDFN
  ENDFN
  /* Check if we meet another patch in the same component */
  IF (NNType EQ rMPG) AND (ComponentId NEQ 0) AND (ComponentId NEQ currPatchId)
     // Find the patch nodes
     nPos2 = PosList[ComponentId]
     // If these are not directly connected
     IF (!LINKED(nnGraph[currType], patchPos, nPos2, 0))
//         d = DISTANCE(NNLoc, startLocation)
        rowDiff = MAX(0, |ROW(NNLoc) - ROW(startLocation)| - 1)
        colDiff = MAX(0, |COL(NNLoc) - COL(startLocation)| - 1)
        d = (rowDiff^2 + colDiff^2)^0.5
        // Create a link
        edgeVar[rEdgeType] = rMPG
        edgeVar[rPatchType] = PatchLayer
        edgeVar[rNode1Id] = currPatchId
        edgeVar[rNode2Id] = ComponentId
        edgeVar[rDist] = d
        Area1 = GET(nnGraph[currType], patchPos, rArea)
        Area2 = GET(nnGraph[currType], nPos2, rArea)
        edgeVar[rWeight] = (Area1 * Area2)^0.5
        edgeVar[rStartLoc] = startLocation
        edgeVar[rEndLoc] = NNLoc
        INSERT LINK(nnGraph[currType], patchPos, nPos2, edgeVar)
        IF (currType EQ 1)
           OVER REGION VECTOR(startLocation,NNLoc)
              PatchLinks = edgeVar[rEdgeType]+1
           ENDFN
        ENDFN
     ENDFN
  ENDFN
  TRANSITIONS = CLASSIFY(NNType) 
                  rNN: (ComponentId EQ 0) AND (PatchesToConnect > 0)
                  rMST: (ComponentId EQ 0) AND (ComponentsToConnect > 0) AND (currPatchId <= NP)
                  rMPG: (ComponentId EQ 0) // Continue until the entire layer has been visited
                ENDFN
  ComponentId = currPatchId
  NNLoc  = startLocation
ENDTR
SPREADTIMESTEP
  IF (startLocation NEQ Location) THEN 0.0000001*(1/MaxDist) ELSE -2
ENDST
/* Patches include diagonal neighbours, but if we spread */
/* using diagonal, then we have problems */
SPREADLOCATION
  IF (PatchId EQ currPatchId)
     maxD = 1.5
  ELSE
     rowDiff = | FLOOR(Location / NUMCOLS) - FLOOR(startLocation / NUMCOLS) | + 1
     colDiff = | (Location % NUMCOLS) - (startLocation % NUMCOLS) | + 1
     maxD = IF ((numSteps + 2 - CEIL((rowDiff^2 + colDiff^2)^0.5)) > 1) THEN 1.5 ELSE 1
  ENDFN
 REGION CENTRED(1, maxD, EUCLIDEAN)
//      DECISION (PatchId NEQ currPatchId) AND (!onEdge OR (PatchLayer EQ currType) OR (DISTANCE(Location, SOURCE     Location) EQ 1))
     DECISION IF (PatchId EQ currPatchId) THEN (ComponentId EQ 0) AND (0 <= PatchLayer < (MinPType + NumPTypes)) AND   (StudyArea > 0)
              ELSE IF (NNType NEQ rMPG) THEN (effectivePatchId[ComponentId] NEQ effectivePatchId[currPatchId]) AND  (0 <= PatchLayer < (MinPType + NumPTypes)) AND (StudyArea > 0)
              ELSE (ComponentId NEQ currPatchId) AND (0 <= PatchLayer < (MinPType + NumPTypes)) AND (StudyArea > 0)
ENDSL
SPREADPROB
  SPREADPROB = 1
//   onEdge = (EdgeInterior EQ Edge) AND (PatchLayer EQ currType)
  IF (PatchId EQ currPatchId)
     startLocation = Location
     numSteps = 0
  ELSE
     startLocation = SOURCE startLocation
     numSteps = SOURCE numSteps + 1
  ENDFN
ENDSP

Suggested Experiments

To explore this cellular automata model further, try the following: