LSStats

Download Model
Download the .scn, .sel and .lse files by clicking on the following link:

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: