### Install sjSDM from CRAN

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

This is the recommended installation method for most users. It installs the package directly from CRAN.

```r
install.packages("sjSDM")
```

--------------------------------

### Install sjSDM Python Package

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/sjSDM/inst/python/README.md

Install the sjSDM Python package using pip. Ensure PyTorch is installed beforehand.

```bash
pip install sjSDM_py
```

--------------------------------

### Interactive Program Startup Notice

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/LICENSE.md

This notice should be displayed when a program starts in interactive mode. It includes copyright information, a statement of no warranty, and details on redistribution rights.

```text
sjSDM Copyright (C) 2020 Maximilian Pichler
This program comes with ABSOLUTELY NO WARRANTY; for details type 'show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type 'show c' for details.
```

--------------------------------

### Install sjSDM Package

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Installs a specific version of the sjSDM package from a local tarball. Ensure the correct version is specified.

```R
install.packages("sjSDM_packages/sjSDM_<version>.tar.gz", repos = NULL, type="source")
```

--------------------------------

### Install Development Version of sjSDM from GitHub

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

Installs the current development version of the sjSDM package directly from the GitHub repository. Dependencies should be installed separately.

```r
devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "master")
```

--------------------------------

### Parallel Processing Setup

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Optimizer_comparison.md

Prepares data for parallel processing and exports necessary variables to cluster nodes. This involves splitting the test data and setting up the parallel environment.

```r
n = 12
cuts = cut(1:nrow(test), breaks = n)
levels(cuts) = 1:n

tests = lapply(1:n, function(i) test[cuts==paste0(i),])

library(snow)
cl = snow::makeCluster(n)
nodes = unlist(snow::clusterEvalQ(cl, paste(Sys.info()[['nodename']], Sys.getpid(), sep='-')))
control = snow::clusterEvalQ(cl, {library(sjSDM)})
snow::clusterExport(cl, list("communities", "test", "tests", "nodes"))
```

--------------------------------

### Install sjSDM Dependencies (GPU or CPU)

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

Automatically installs the necessary dependencies for the sjSDM package. Specify 'gpu' or 'cpu' for the desired version.

```r
sjSDM::install_sjSDM(version = "gpu") # or
sjSDM::install_sjSDM(version = "cpu")
```

--------------------------------

### Load s-jsdm and Set Seed

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Internal.md

Loads the s-jsdm library and sets the random seed for reproducibility. This is a common starting point for simulations.

```r
library(sjSDM)
set.seed(42)
```

--------------------------------

### Build and Train a Linear JSDM Model

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/sjSDM/inst/python/README.md

Example demonstrating how to build, train, and inspect a linear Joint Species Distribution Model using sjSDM. Requires PyTorch and NumPy. Consider setting 'df' to 'number of species / 2' and use 'bias=True' in Layer_dense for species intercepts.

```python
import sjSDM_py as sa
import numpy as np
Env = np.random.randn(100, 5)
Occ = np.random.binomial(1, 0.5, [100, 10])

model = sa.Model_base(5) # input_shape == number of environmental predictors
model.add_layer(sa.layers.Layer_dense(hidden=10)) # number of hidden units in the layer == number of species
model.build(df=5, optimizer=sa.optimizer_adamax(lr=0.1, weight_decay = 0.01)) # df = degree of freedom 
model.fit(X = Env, Y = Occ)
print(model.weights_numpy)
print(model.get_cov())
```

--------------------------------

### Load s-jsdm Package

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

Loads the s-jsdm package into the R session. Ensure the package is installed before running this command.

```r
library(sjSDM)
```

--------------------------------

### Importing Libraries for sjSDM Simulation

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/EM sjSDM.ipynb

Imports necessary libraries for simulation, including PyTorch, NumPy, Matplotlib, Statsmodels, and SciPy. Ensure these libraries are installed before running the simulation.

```python
import torch
import numpy as np
import numpy as np, numpy.linalg
import matplotlib.pyplot as plt
from dataclasses import dataclass
import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
from scipy.optimize import minimize
from dataclasses import dataclass
import numpy as np
```

--------------------------------

### Run Process-Based Model Simulation Script

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Execute the R script for simulating data from the process-based model. Ensure sjSDM version 0.1.8 is installed.

```r
source("analysis/scripts/7_simulate_from_process_based_model.R")
```

--------------------------------

### Install PyTorch for sjSDM 0.0.7.9000

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Installs compatible PyTorch versions (1.8.0) and related libraries using Conda for sjSDM version 0.0.7.9000. This is a manual process to avoid conflicts with newer PyTorch versions.

```bash
conda create -n r-reticulate python=3.6
conda activate r-reticulate
conda install -n r-reticulate  pytorch==1.8.0 torchvision==0.9.0 torchaudio==0.8.0 cudatoolkit=10.2 -c pytorch
python -m pip install tqdm torch_optimizer pyro-ppl==1.6.0
```

--------------------------------

### Fit Model with CAR=False

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/EM sjSDM.ipynb

Calls the `fit_model` function with the `CAR` parameter set to `False`. This is a typical usage example for fitting the model without the CAR component.

```python
fit_model(data, CAR=False, det=True)
```

--------------------------------

### Fit MVP jSDM Model

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/EM sjSDM.ipynb

Fits a multivariate probit (MVP) jSDM using PyTorch. This function handles model setup, optimization, and likelihood calculation for both standard and CAR random effects.

```python
def fit_model(data, 
              outer_epochs=12, 
              inner_epochs=3, 
              device = "cpu:0", 
              optim = "torch", 
              likelihood_type="mvp", 
              det=True,
              CAR=False
             ):
    E = data.X.shape[1]
    N = data.X.shape[0]
    SP = data.Y.shape[1]
    X, Y, G, indices, D = data.X, data.Y, data.G, data.g, data.D
    if CAR is True:
        G = G**2
    dev = torch.device(device)
    XT = torch.tensor(X, dtype=torch.float32, device=dev)
    YT = torch.tensor(Y, dtype=torch.float32, device=dev)
    torch.autograd.set_detect_anomaly(True)
    W = torch.tensor(np.random.normal(0., 0.001, size=(XT.shape[1],Y.shape[1])), dtype=torch.float32, requires_grad=True, device=dev)
    r_dim = Y.shape[1]
    df = int(np.rint(Y.shape[1]/2))
    low = -np.sqrt(6.0/(r_dim+df)) # type: ignore
    high = np.sqrt(6.0/(r_dim+df)) # type: ignore      
    sigma = torch.tensor(np.random.uniform(low, high, [r_dim, df]), requires_grad = True, dtype=torch.float32, device=dev) # type: ignore

    @torch.jit.script
    def likelihood(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
        noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)
        E = torch.sigmoid(   torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu).mul(alpha)   ).mul(0.999999).add(0.0000005)
        logprob = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys)).neg().sum(dim = 2).neg()
        maxlogprob = logprob.max(dim = 0).values
        Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)
        return Eprob.log().neg().sub(maxlogprob)

    scale_log = torch.tensor(np.random.normal(0.0,0.001, [1]), dtype=torch.float32, requires_grad=True, device=dev)
    res = torch.tensor(np.random.normal(0.0,0.001, [G, 1]), dtype=torch.float32, requires_grad=True, device=dev)
    zero_intercept = torch.zeros([1], dtype=torch.float32, device=dev)
    zero_CAR = torch.zeros([G], dtype=torch.float32, device=dev)

    opt1 = torch.optim.RMSprop([W, scale_log, sigma],lr=0.01)

    if CAR is True:
        D = torch.tensor(D, dtype=torch.float32, device=dev)
    opt2 = torch.optim.LBFGS([res], lr = 0.1)
    
    const_val = torch.tensor(0.0, dtype=torch.float32, device=dev)
    const_cov = torch.eye(G, dtype=torch.float32, device=dev)*0.01

    soft = lambda t: torch.nn.functional.softplus(t)+0.0001

    indices_T = torch.tensor(indices, dtype=torch.long)
    
    def re_loss():
        if CAR is not True:
            return -torch.distributions.Normal(zero_intercept, soft(scale_log)).log_prob(res).sum()
        else:
            return -torch.distributions.MultivariateNormal(zero_CAR, (-soft(scale_log)*D).exp()+const_cov).log_prob(res.reshape([1,-1])).sum()

    
    def ll(res, W,sigma, XT, YT, indices_T):
        pred = XT@W+res[indices_T,:]#*scale_log.exp()
        #loss = -torch.distributions.Normal(loc=pred, scale=soft(scale_log_2) ).log_prob(YT).sum()
        if likelihood_type == "mvp":
            loss = likelihood(pred, YT, sigma, XT.shape[0], 100, df, 1.7012, device, torch.float32).sum()
        else:
            loss = -torch.distributions.Binomial(total_count=1, probs= torch.sigmoid(pred*1.7012) ).log_prob(YT).sum()

```

--------------------------------

### Standard Copyright Notice for Source Files

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/LICENSE.md

This is the recommended copyright notice to be placed at the beginning of each source file. It includes the program name, copyright holder, licensing terms (GNU GPL v3 or later), and warranty disclaimer.

```text
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) 2020 Maximilian Pichler

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
```

--------------------------------

### Get Citation Information for sjSDM

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

Type this command in R to get citation information for the sjSDM package when using it for research.

```r
citation("sjSDM")
```

--------------------------------

### Initialize and Build sjSDM Model in Python

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/README.md

Demonstrates how to initialize, add environmental data, and build an sjSDM model in Python using the sjSDM_py package. It specifies the device, data types, optimizer, and scheduler.

```python
import sjSDM_py as fa
import numpy as np
import torch

Env = np.random.randn(100, 5)
Occ = np.random.binomial(1, 0.5, [100, 10])

model = fa.Model_sjSDM(device=torch.device("cpu"), dtype=torch.float32)
model.add_env(5, 10)
model.build(5, optimizer=fa.optimizer_adamax(0.001),scheduler=False)
model.fit(Env, Occ, batch_size = 20, epochs = 10)
```

--------------------------------

### Define Optimization Parameters

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Optimizer_comparison.md

Sets up a comprehensive grid of parameters for optimizer comparison, including optimizer type, learning rate, epochs, batch size, and number of species.

```r
optims = c("adamax", "rmsprop", "accsgd", "adabound", "sgd", "diffgrad")
scheduler = c(TRUE, FALSE)
lr = c(0.01, 0.005, 0.002, 0.001)
epochs = c(100, 200)
bs = c(5, 10, 30)
sp = c(25, 50, 100)
iter = 1:10

test = data.frame(expand.grid(optims, scheduler, lr, epochs, bs, sp, iter))
colnames(test) = c("optim", "scheduler","lr","epochs", "bs", "sp", "iter")
test$acc = NA
test$grads = NA
test$rmse = NA
test = test[order(test$sp),]
```

--------------------------------

### Clone TestingHMSC Repository

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Clone the testingHMSC GitHub repository to your working directory. This is a prerequisite for running the process-based simulations.

```bash
git clone https://github.com/javirudolph/testingHMSC
```

--------------------------------

### Load Required Libraries

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

Loads the necessary libraries for running the simulation experiments, including sjSDM and Hmsc.

```r
set.seed(42)
library(sjSDM)
library(Hmsc)
library(snow)
```

--------------------------------

### Benchmark GLLVM Package

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Runs the GLLVM package for performance benchmarking, comparing its runtime against sjSDM.

```R
source("analysis/scripts/1_gllvm.R)"
```

--------------------------------

### Accessing data.W in LA.ipynb

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA.ipynb

Demonstrates how to access the 'W' attribute from a 'data' object, likely retrieving a weight matrix or similar data structure.

```python
data.W
```

--------------------------------

### Simulate Species Data with Spatial Effects

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Internal.md

Generates simulated species occurrence data (Y) influenced by environmental factors (Env), spatial coordinates (XY), and a species-specific covariance matrix (Sigma). This setup is used to test different spatial modeling techniques.

```r
SP = 10
Sigma = diag(1.0, SP)
Sigma[] = 0.9
Sigma[1:5,] = Sigma[,1:5] = 0.0
diag(Sigma) = 1.0
fields::image.plot(Sigma)
beta = c(rep(1.5, 5), rep(0.3, 5))
Env = rnorm(500)
XY = matrix(rnorm(1000), 500, 2)
betaSP = rep(0.3, SP)
Y = 1*((Env %*% t(beta) + (XY[,1,drop=FALSE]*XY[,2,drop=FALSE]) %*% t(betaSP) +  mvtnorm::rmvnorm(500, sigma = Sigma))>0)
```

--------------------------------

### Run s-jsdm and Hmsc simulations for 25 species

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

This snippet sets up and runs simulations for 25 species, comparing s-jsdm and Hmsc models. It generates data, fits models, and calculates performance metrics like correlation accuracy and RMSE.

```r
data = lapply(1:5, function(l) create(env = 3L, n = 500L, sp = 25L, l = l))
results_3 = matrix(NA, 5, 5) # for acc/rmse/acc/rmse/n_latent
for(i in 1:5) {
  sjSDM = sjSDM(data[[i]]$X, data[[i]]$Y, formula = ~0+., step_size = 20L, 
                iter = 100L, device = 0L, sampling = 500L, learning_rate = 0.005)
  sp_sjSDM = getCov(sjSDM)
  results_3[i,1] = data[[i]]$corr_acc(sp_sjSDM)
  results_3[i,2] = sqrt(mean(as.vector(coef(sjSDM)[[1]] - data[[i]]$SPW)^2))
}


cl = snow::makeCluster(2L)
snow::clusterExport(cl, list("data"))
ev = snow::clusterEvalQ(cl, {
  library(Hmsc) 
  })

hmsc_res = 
  snow::parLapply(cl, 1:5, function(i) {
  
    hmsc = list()
    studyDesign = data.frame(sample = as.factor(1:nrow(data[[i]]$Y)))
    rL = HmscRandomLevel(units = studyDesign$sample)
    model = Hmsc(Y = data[[i]]$Y, XData = data.frame(data[[i]]$X), XFormula = ~0 + .,
                 studyDesign = studyDesign, ranLevels = list(sample = rL), distr = "probit")
    model = sampleMcmc(model, thin = 50, samples = 500, transient = 5000,verbose = 5000,
                       nChains = 1L) # 50,000 iterations
    correlation = computeAssociations(model)[[1]]$mean
    beta = Hmsc::getPostEstimate(model, "Beta")$mean
    return(c(data[[i]]$corr_acc(correlation),
             sqrt(mean(as.vector(beta - data[[i]]$SPW)^2)) ,
             dim(Hmsc::getPostEstimate(model, "Lambda")$mean)[1]))

})
snow::stopCluster(cl)
results_3[,3:5] = abind::abind(hmsc_res, along = 0L)
```

--------------------------------

### Fit JSDM Models (Process-Based)

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Run R scripts to fit JSDM models to the simulated data from the process-based model. Two scripts are provided for different fitting configurations (A and B).

```r
source("analysis/scripts/7_process_based_models_A.R")
```

```r
source("analysis/scripts/7_process_based_models_B.R")
```

--------------------------------

### PyTorch Model Training Loop

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA.ipynb

This snippet shows a training loop for a PyTorch model, including loss calculation, backpropagation, and optimizer steps. It also includes conditional printing of loss values and returns model parameters.

```python
D_tmp = hess.index_select(0, ind2).index_select(1, ind2)
const_val = torch.eye(ind2.shape[0], device=dev, dtype=torch.float32)*0.01
logDA=(D_tmp+const_val).inverse().logdet()
loss2 = ((ind2.shape[0]*0.5*torch.log((2*torch.tensor(3.14, dtype=torch.float32))) - 0.5*logDA))/adapt/x.shape[0]
loss+=loss2
loss = loss
loss.backward()
optimizer.step()

optimizer.zero_grad()



if i % 20 == 0:
     print([loss.item(), loss2])
     
return [(torch.exp(scale_log)).cpu().data.numpy().tolist(), 
            (torch.exp(scale_log_normal)).cpu().data.numpy().tolist(), 
            W.cpu().data.numpy()]
```

--------------------------------

### Prepare eDNA Data

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Prepares the eDNA Fungi dataset for analysis. Note: the path argument within the script may need to be adjusted.

```R
source("data/eDNA/prepare_eDNA.R") # please change the path
```

--------------------------------

### Run Case Study 1 with sjSDM

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Applies sjSDM to the datasets from Wilkinson et al. 2019, likely for a specific case study analysis.

```R
source("analysis/scripts/4_case_study_1.R")
```

--------------------------------

### Simulate Communities

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Optimizer_comparison.md

Generates simulated ecological community data with varying numbers of species.

```r
communities = lapply(c(25, 50, 100), function(sp) lapply(1:10, function(i) simulate_SDM(env = 3L, sites = 100L, species = sp)))
```

--------------------------------

### Tensor Outputs from LA.ipynb

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA.ipynb

Displays a series of tensor values with gradient function information, indicating a computation or optimization process.

```text
tensor(186.3492, grad_fn=<LogdetBackward0>)
tensor(186.2140, grad_fn=<LogdetBackward0>)
tensor(186.0867, grad_fn=<LogdetBackward0>)
tensor(185.9677, grad_fn=<LogdetBackward0>)
tensor(185.8521, grad_fn=<LogdetBackward0>)
tensor(185.7405, grad_fn=<LogdetBackward0>)
tensor(185.6338, grad_fn=<LogdetBackward0>)
tensor(185.5270, grad_fn=<LogdetBackward0>)
tensor(185.4211, grad_fn=<LogdetBackward0>)
tensor(185.3172, grad_fn=<LogdetBackward0>)
tensor(185.2108, grad_fn=<LogdetBackward0>)
tensor(185.1030, grad_fn=<LogdetBackward0>)
tensor(184.9951, grad_fn=<LogdetBackward0>)
tensor(184.8828, grad_fn=<LogdetBackward0>)
tensor(184.7678, grad_fn=<LogdetBackward0>)
tensor(184.6514, grad_fn=<LogdetBackward0>)
tensor(184.5295, grad_fn=<LogdetBackward0>)
tensor(184.4041, grad_fn=<LogdetBackward0>)
tensor(184.2766, grad_fn=<LogdetBackward0>)
tensor(184.1432, grad_fn=<LogdetBackward0>)
tensor(184.0059, grad_fn=<LogdetBackward0>)
tensor(183.8665, grad_fn=<LogdetBackward0>)
tensor(183.7212, grad_fn=<LogdetBackward0>)
tensor(183.5721, grad_fn=<LogdetBackward0>)
tensor(183.4210, grad_fn=<LogdetBackward0>)
tensor(183.2644, grad_fn=<LogdetBackward0>)
tensor(183.1045, grad_fn=<LogdetBackward0>)
tensor(182.9430, grad_fn=<LogdetBackward0>)
tensor(182.7766, grad_fn=<LogdetBackward0>)
tensor(182.6073, grad_fn=<LogdetBackward0>)
tensor(182.4371, grad_fn=<LogdetBackward0>)
tensor(182.2626, grad_fn=<LogdetBackward0>)
tensor(182.0859, grad_fn=<LogdetBackward0>)
tensor(181.9089, grad_fn=<LogdetBackward0>)
tensor(181.7283, grad_fn=<LogdetBackward0>)
tensor(181.5461, grad_fn=<LogdetBackward0>)
tensor(181.3644, grad_fn=<LogdetBackward0>)
tensor(181.1796, grad_fn=<LogdetBackward0>)
tensor(180.9939, grad_fn=<LogdetBackward0>)
tensor(180.8093, grad_fn=<LogdetBackward0>)
tensor(180.6221, grad_fn=<LogdetBackward0>)
tensor(180.4346, grad_fn=<LogdetBackward0>)
tensor(180.2487, grad_fn=<LogdetBackward0>)
tensor(180.0608, grad_fn=<LogdetBackward0>)
tensor(179.8730, grad_fn=<LogdetBackward0>)
tensor(179.6872, grad_fn=<LogdetBackward0>)
tensor(179.4998, grad_fn=<LogdetBackward0>)
tensor(179.3128, grad_fn=<LogdetBackward0>)
tensor(179.1282, grad_fn=<LogdetBackward0>)
tensor(178.9422, grad_fn=<LogdetBackward0>)
tensor(178.7569, grad_fn=<LogdetBackward0>)
tensor(178.5742, grad_fn=<LogdetBackward0>)
tensor(178.3903, grad_fn=<LogdetBackward0>)
tensor(178.2073, grad_fn=<LogdetBackward0>)
tensor(178.0269, grad_fn=<LogdetBackward0>)
tensor(177.8454, grad_fn=<LogdetBackward0>)
tensor(177.6649, grad_fn=<LogdetBackward0>)
```

--------------------------------

### Benchmark BayesComm Package

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Runs the BayesComm package for performance benchmarking, comparing its runtime against sjSDM.

```R
source("analysis/scripts/1_BayesCommDiag.R")
```

--------------------------------

### Simulate s-jsdm Data

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/EM sjSDM.ipynb

Simulates data for the s-jsdm model with specified parameters for species (SP) and number of observations (N).

```python
data = simulate(SP = 100, N = 1000  )
```

--------------------------------

### Run s-jsdm and Hmsc simulations for 50 species

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

This snippet executes simulations for 50 species, comparing s-jsdm and Hmsc models. It generates data, fits models, and calculates performance metrics such as correlation accuracy and RMSE.

```r
data = lapply(1:5, function(l) create(env = 3L, n = 500L, sp = 50L, l = l))
results_4 = matrix(NA, 5, 5) # for acc/rmse/acc/rmse/n_latent
for(i in 1:5) {
  sjSDM = sjSDM(data[[i]]$X, data[[i]]$Y, formula = ~0+., step_size = 20L, 
                iter = 100L, device = 0L, sampling = 500L, learning_rate = 0.005)
  sp_sjSDM = getCov(sjSDM)
  results_4[i,1] = data[[i]]$corr_acc(sp_sjSDM)
  results_4[i,2] = sqrt(mean(as.vector(coef(sjSDM)[[1]] - data[[i]]$SPW)^2))
}


cl = snow::makeCluster(2L)
snow::clusterExport(cl, list("data"))
ev = snow::clusterEvalQ(cl, {
  library(Hmsc) 
  })

hmsc_res = 
  snow::parLapply(cl, 1:5, function(i) {
  
    hmsc = list()
    studyDesign = data.frame(sample = as.factor(1:nrow(data[[i]]$Y)))
    rL = HmscRandomLevel(units = studyDesign$sample)
    model = Hmsc(Y = data[[i]]$Y, XData = data.frame(data[[i]]$X), XFormula = ~0 + .,
                 studyDesign = studyDesign, ranLevels = list(sample = rL), distr = "probit")
    model = sampleMcmc(model, thin = 50, samples = 500, transient = 5000,verbose = 5000,
                       nChains = 1L) # 50,000 iterations
    correlation = computeAssociations(model)[[1]]$mean
    beta = Hmsc::getPostEstimate(model, "Beta")$mean
    return(c(data[[i]]$corr_acc(correlation),
             sqrt(mean(as.vector(beta - data[[i]]$SPW)^2)) ,
             dim(Hmsc::getPostEstimate(model, "Lambda")$mean)[1]))

})
snow::stopCluster(cl)
results_4[,3:5] = abind::abind(hmsc_res, along = 0L)
```

--------------------------------

### Load Libraries

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Optimizer_comparison.md

Loads necessary libraries for sjSDM and data manipulation.

```r
set.seed(42)
library(sjSDM)
library(dplyr)
```

--------------------------------

### Model Fitting with s-jsdm

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA2.ipynb

Fits a statistical model using the s-jsdm library. Specify the data, device for computation, whether to compute the determinant, batch size, and number of epochs for training.

```python
fit_model(data, device="cpu:0", det=True, batch_size=10, epochs=100)
```

--------------------------------

### Data Simulation and MixedLM Fitting with Statsmodels

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA.ipynb

This snippet simulates data for a mixed-effects model and then fits the model using Statsmodels' `mixedlm` function. It prints the summary of the fitted model.

```python
data = simulate(N=1000, G = 50, sd_re=2.3)
print(2.3**2)
dataset = pd.DataFrame(np.concatenate([data.X, data.Y, np.reshape(data.g, [data.g.shape[0],1])], axis = 1),
                       columns = ["X1", "X2", "Y", "g"])
md = smf.mixedlm("Y~0+X1+X2", dataset, groups=dataset["g"], )
mdf = md.fit(reml=False)
print(mdf.summary())
```

--------------------------------

### Run Sparse GLLVM

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Executes the GLLVM package for scenarios involving sparse species-species associations.

```R
source("analysis/scripts/6_sparse_gllvm.R)"
```

--------------------------------

### Result List from LA.ipynb

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA.ipynb

Presents a list containing numerical values and a NumPy array, likely representing intermediate or final results of a calculation.

```text
[1.9286086559295654,
 1.0407851934432983,
 array([[-1.3586684 ],
        [ 0.22159573]], dtype=float32)]
```

--------------------------------

### Run s-jsdm and Hmsc simulations for 100 species

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

This snippet conducts simulations for 100 species, comparing s-jsdm and Hmsc models. It generates data, fits models, and calculates performance metrics including correlation accuracy and RMSE.

```r
data = lapply(1:5, function(l) create(env = 3L, n = 500L, sp = 100L, l = l))
results_5 = matrix(NA, 5, 5) # for acc/rmse/acc/rmse/n_latent
for(i in 1:5) {
  sjSDM = sjSDM(data[[i]]$X, data[[i]]$Y, formula = ~0+., step_size = 20L, 
                iter = 100L, device = 0L, sampling = 500L, learning_rate = 0.005)
  sp_sjSDM = getCov(sjSDM)
  results_5[i,1] = data[[i]]$corr_acc(sp_sjSDM)
  results_5[i,2] = sqrt(mean(as.vector(coef(sjSDM)[[1]] - data[[i]]$SPW)^2))
}


cl = snow::makeCluster(2L)
snow::clusterExport(cl, list("data"))
ev = snow::clusterEvalQ(cl, {
  library(Hmsc) 
  })

hmsc_res = 
  snow::parLapply(cl, 1:5, function(i) {
  
    hmsc = list()
    studyDesign = data.frame(sample = as.factor(1:nrow(data[[i]]$Y)))
    rL = HmscRandomLevel(units = studyDesign$sample)
    model = Hmsc(Y = data[[i]]$Y, XData = data.frame(data[[i]]$X), XFormula = ~0 + .,
                 studyDesign = studyDesign, ranLevels = list(sample = rL), distr = "probit")
    model = sampleMcmc(model, thin = 50, samples = 500, transient = 5000,verbose = 5000,
                       nChains = 1L) # 50,000 iterations
    correlation = computeAssociations(model)[[1]]$mean
    beta = Hmsc::getPostEstimate(model, "Beta")$mean
    return(c(data[[i]]$corr_acc(correlation),
             sqrt(mean(as.vector(beta - data[[i]]$SPW)^2)) ,
             dim(Hmsc::getPostEstimate(model, "Lambda")$mean)[1]))

})
snow::stopCluster(cl)
results_5[,3:5] = abind::abind(hmsc_res, along = 0L)
```

--------------------------------

### Generate Data for Analysis

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Executes an R script to generate datasets required for various analyses, including runtime benchmarking.

```R
source("./analysis/scripts/1_generate_data.R")
```

--------------------------------

### Benchmark sjSDM on CPU

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Runs the sjSDM model on the CPU for performance benchmarking, specifically for dense association scenarios.

```R
source("analysis/scripts/1_cpu_sjSDM.R")
```

--------------------------------

### Simulation Test 1: 5 Species, 3 Env Covariates, 500 Sites

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

Runs a simulation with 5 species and 3 environmental covariates across 500 sites, testing HmscLVM with varying numbers of latent variables. It compares sjSDM and Hmsc models.

```r
data = lapply(1:5, function(l) create(env = 3L, n = 500L, sp = 5L, l = l))
results_1 = matrix(NA, 5, 5) # for acc/rmse/acc/rmse/n_latent
for(i in 1:5) {
  sjSDM = sjSDM(data[[i]]$X, data[[i]]$Y, formula = ~0+., step_size = 20L, 
                iter = 100L, device = 0L, sampling = 500L, learning_rate = 0.005)
  sp_sjSDM = getCov(sjSDM)
  results_1[i,1] = data[[i]]$corr_acc(sp_sjSDM)
  results_1[i,2] = sqrt(mean(as.vector(coef(sjSDM)[[1]] - data[[i]]$SPW)^2))
}


cl = snow::makeCluster(2L)
snow::clusterExport(cl, list("data"))
ev = snow::clusterEvalQ(cl, {
  library(Hmsc)
  })

hmsc_res = 
  snow::parLapply(cl, 1:5, function(i) {
  
    hmsc = list()
    studyDesign = data.frame(sample = as.factor(1:nrow(data[[i]]$Y)))
    rL = HmscRandomLevel(units = studyDesign$sample)
    model = Hmsc(Y = data[[i]]$Y, XData = data.frame(data[[i]]$X), XFormula = ~0 + .,
                 studyDesign = studyDesign, ranLevels = list(sample = rL), distr = "probit")
    model = sampleMcmc(model, thin = 50, samples = 500, transient = 5000,verbose = 5000,
                       nChains = 1L) # 50,000 iterations
    correlation = computeAssociations(model)[[1]]$mean
    beta = Hmsc::getPostEstimate(model, "Beta")$mean
    return(c(data[[i]]$corr_acc(correlation),
             sqrt(mean(as.vector(beta - data[[i]]$SPW)^2)) ,
             dim(Hmsc::getPostEstimate(model, "Lambda")$mean)[1]))

})

snow::stopCluster(cl)
results_1[,3:5] = abind::abind(hmsc_res, along = 0L)
```

--------------------------------

### Simulate from LVM and Fit sjSDM

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Execute the R script to simulate data from a Latent Variable Model (LVM) and then fit sjSDM, Hmsc, and gllvm models to this simulated data. Requires sjSDM version 0.1.8.

```r
source("analysis/scripts/8_LVM_simulation.R")
```

--------------------------------

### PyTorch Mixed-Effects Model Fitting

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/LA2.ipynb

Fits a mixed-effects model using PyTorch, optimizing parameters like regression coefficients, scale, and random effects. Requires data in Simulation format.

```python
#@ray.remote(num_gpus=1,num_cpus=2)
def fit_model(data: Simulation, det=True, CAR=False, device = "cuda:0", batch_size = 100, epochs = 40) -> list:
     N, E = data.X.shape
     SP = data.Y.shape[1]
     X, Y, G, indices = data.X, data.Y, data.G, data.g
     dev = torch.device(device)
    
     XT = torch.tensor(X, dtype=torch.float32, device=torch.device("cpu:0"))
     YT = torch.tensor(Y, dtype=torch.float32, device=torch.device("cpu:0"))
     indices_T = torch.tensor(indices, dtype=torch.long, device=torch.device("cpu:0"))

     # Variables
     W = torch.tensor(np.random.normal(0.0,0.001, [XT.shape[1], YT.shape[1]]), dtype=torch.float32, device=dev, requires_grad=True)
     scale_log = torch.tensor(1.0, dtype=torch.float32, requires_grad=True, device=dev)
     scale_log_normal = torch.tensor(1.0, dtype=torch.float32, requires_grad=True, device=dev)
     res = torch.tensor(np.random.normal(0.0,0.001, [G, 1]), dtype=torch.float32, requires_grad=True, device=dev)
     soft = lambda t: torch.nn.functional.softplus(t)+0.0001
     zero_intercept = torch.zeros([1], dtype=torch.float32, device=dev)
     loss2 = torch.zeros([1], dtype=torch.float32, device=dev)

     
     adapt = torch.tensor(np.rint(XT.shape[0]/batch_size).tolist(), dtype=torch.float32, device=dev)
     const = torch.tensor(np.log(2*(np.pi)), dtype=torch.float32, device=dev)
     
     def ll(res, W, XT, YT, indices_T, scale_log, scale_log_normal):
        pred = XT@W+res[indices_T,:]
        loss = -torch.distributions.Normal(pred, torch.clamp(scale_log_normal, 0.00001, 100.)).log_prob(YT).sum()/XT.shape[0]
        loss += -torch.distributions.Normal(zero_intercept, torch.clamp(scale_log,0.0001, 100.)).log_prob(res[indices_T.unique()]).sum()/adapt/XT.shape[0]#/G/indices_T.unique().shape[0] #/XT.shape[0]
        #loss += ((res[indices_T.unique()].pow(2.0))*(sigma_res)*0.5).sum()/res[indices_T.unique()].shape[0]#/crit_factor
        return loss


     optimizer = torch.optim.Adamax([W, scale_log,scale_log_normal, res], lr = 0.1)

     dataset = torch.utils.data.TensorDataset(XT, YT, indices_T)
     dataLoader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True)
     

     for i in range(epochs):

        for x, y, inds in dataLoader:
             optimizer.zero_grad()
             loss = ll(res, W, x.to(dev), y.to(dev), inds.to(dev),scale_log, scale_log_normal)#/x.shape[0]
             if det is True:

                  hess = torch.autograd.functional.hessian(lambda res: ll(res, W, x.to(dev), y.to(dev), inds.to(dev),scale_log, scale_log_normal), res, create_graph=True).squeeze()
                  ind2 = inds.unique()
                  
                  D_tmp = hess.index_select(0, ind2).index_select(1, ind2)
                  const_val = torch.eye(ind2.shape[0], device=dev, dtype=torch.float32)*0.01
                  logDA=(D_tmp+const_val).inverse().logdet()/G/indices_T.unique().shape[0]
                  loss2 = -(0.5*logDA)/adapt/x.shape[0]
                  loss+=loss2 + const*G/2.
             loss = loss
             loss.backward()
             optimizer.step()          


        if i % 2 == 0:
            print([loss.item(), loss2, scale_log.item()])
               
     return [(scale_log).cpu().data.numpy().tolist(), 
             (scale_log_normal).cpu().data.numpy().tolist(), 
             W.cpu().data.numpy()]

```

--------------------------------

### Run Optimizer Comparison in Parallel

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/Optimizer_comparison.md

Executes the sjSDM model training with different optimizers and parameters in parallel across multiple cluster nodes. It calculates accuracy, gradient norms, and RMSE for each configuration.

```r
results = parLapply(cl, 1:n, function(n) {
  tmp = tests[[n]]
  for(i in 1:nrow(tmp)) {
    if(tmp[i,]$sp == 25) j = 1
    else if(tmp[i,]$sp == 50) j = 2
    else j = 3
    com = communities[[j]][[tmp[i,]$iter]]
    
    opt = 
    switch (as.character(tmp[i,]$optim), 
      adamax = Adamax(),
      rmsprop = RMSprop(),
      accsgd = AccSGD(),
      adabound = AdaBound(),
      sgd = SGD(),
      diffgrad = DiffGrad()
    )
    
    if(n %in% 1:4) device = 0
    else if(n %in% 5:8) device = 1
    else if(n %in% 9:12) device = 2
    
    
    m = sjSDM(com$response, com$env_weights, 
              learning_rate = tmp[i,]$lr, 
              step_size = tmp[i,]$bs, 
              iter = tmp[i,]$epochs, 
              control = sjSDMControl(opt, tmp[i,]$scheduler),
              device = device)
    
    tmp[i,]$acc = com$corr_acc(getCov(m))
    tmp[i,]$grads = m$model$params[[2]][[1]]$grad$pow(2.0)$mean()$sqrt()$item()
    tmp[i,]$rmse = sqrt(mean( ( t(coef(m)[[1]])[-1,] - com$species_weights )^2 ))
    rm(m)
  }
  return(tmp)
  
})
```

--------------------------------

### Run Fungi eDNA Models

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Executes the main analysis models on the prepared Fungi eDNA dataset after tuning.

```R
source("analysis/scripts/5_Fungi_models.R")
```

--------------------------------

### Run Sparse BayesComm

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Executes the BayesComm package for scenarios involving sparse species-species associations.

```R
source("analysis/scripts/6_sparse_bc.R")
```

--------------------------------

### Analyze Covariance Behavior

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/publications/Pichler and Hartig- 2020/README.md

Runs a script to specifically investigate the behavior of covariance estimation in the models.

```R
source("analysis/scripts/3_covariance_behaviour.R")
```

--------------------------------

### Fit s-jsdm Model

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/sjSDM with random effects/EM sjSDM.ipynb

Fits the s-jsdm model using simulated data with specified likelihood function ('binomial'), intercept setting ('species'), degrees of freedom (df), batch size, and determinant calculation flag.

```python
fit_model(data, ll="binomial", intercept="species", df = 40, batch_size = 200, det = False)
```

--------------------------------

### Simulation Test 2: 10 Species, 3 Env Covariates, 500 Sites

Source: https://github.com/theoreticalecology/s-jsdm/blob/master/Code/SimulationExperiments/HmscLVMsimulation.md

Runs a simulation with 10 species and 3 environmental covariates across 500 sites, testing HmscLVM with varying numbers of latent variables. It compares sjSDM and Hmsc models.

```r
data = lapply(1:5, function(l) create(env = 3L, n = 500L, sp = 10L, l = l))
results_2 = matrix(NA, 5, 5) # for acc/rmse/acc/rmse/n_latent
for(i in 1:5) {
  sjSDM = sjSDM(data[[i]]$X, data[[i]]$Y, formula = ~0+., step_size = 20L, 
                iter = 100L, device = 0L, sampling = 500L, learning_rate = 0.005)
  sp_sjSDM = getCov(sjSDM)
  results_2[i,1] = data[[i]]$corr_acc(sp_sjSDM)
  results_2[i,2] = sqrt(mean(as.vector(coef(sjSDM)[[1]] - data[[i]]$SPW)^2))
}

cl = snow::makeCluster(2L)
snow::clusterExport(cl, list("data"))
ev = snow::clusterEvalQ(cl, {
  library(Hmsc)
  })
hmsc_res = 
  snow::parLapply(cl, 1:5, function(i) {
  
    hmsc = list()
    studyDesign = data.frame(sample = as.factor(1:nrow(data[[i]]$Y)))
    rL = HmscRandomLevel(units = studyDesign$sample)
    model = Hmsc(Y = data[[i]]$Y, XData = data.frame(data[[i]]$X), XFormula = ~0 + .,
                 studyDesign = studyDesign, ranLevels = list(sample = rL), distr = "probit")
    model = sampleMcmc(model, thin = 50, samples = 500, transient = 5000,verbose = 5000,
                       nChains = 1L) # 50,000 iterations
    correlation = computeAssociations(model)[[1]]$mean
    beta = Hmsc::getPostEstimate(model, "Beta")$mean
    return(c(data[[i]]$corr_acc(correlation),
             sqrt(mean(as.vector(beta - data[[i]]$SPW)^2)),
             dim(Hmsc::getPostEstimate(model, "Lambda")$mean)[1]))

})
snow::stopCluster(cl)
results_2[,3:5] = abind::abind(hmsc_res, along = 0L)
```