Python example array job

import numpy as np
#import matplotlib.pyplot as plt
import os

##  'Given function, $f(x)=\frac{4}{1+x^2}$ and a density function that is similar to $f(x)$
##  is $p(x)=\\frac{1}{3}(4-2x)$, 
##  remember a density function must be non-negative and $\int_{-\infty}^{\infty} p(x)dx=1$'
##  'The goal is evaluating the definite integral $\\int_0^1 f(x)dx$')

##============================================================================================================
## get the $SLURM_ARRAY_TASK_ID environment variable value
## this value will be utilized to carry out the integration
## with different number of samples to check if number of
## samples makes any difference
##============================================================================================================
if ("SLURM_ARRAY_TASK_ID" in os.environ):
    indx = int(os.environ["SLURM_ARRAY_TASK_ID"])
    print("SLURM_ARRAY_TASK_ID {0}".format(indx))
else:
    print("${SLURM_ARRAY_TASK_ID} is not found, quitting!")
    quit()

x = np.linspace(0,1,20)
fx = lambda x: 4.0/(1+x**2)
px = lambda x: 1/3.0*(4-2*x)
sumv = 0.0;
numsample = 10 * 2**(indx+1)
seed = 200 * (indx+1)
np.random.seed(seed)
sumsv = 0.0

### using variance reduction technique
for ii in range(numsample):
    y = np.random.random()
    sumv += fx(y)/px(y)
    sumsv += (fx(y)/px(y))**2 
integral_val = sumv/numsample
print("Value of the definite integral, assuming a distribution close to original" + \
     " function: {0:2.4f}".format(integral_val))
error_sigma = np.sqrt((sumsv/numsample-(sumv/numsample)**2)) / np.sqrt(numsample)
print("Estimated stdandard deviation: {0:2.6f}".format(error_sigma))
print("Within 3 std: {0:2.3f}-{1:2.3f}". format(integral_val-3.0*error_sigma, integral_val+3.0*error_sigma))
print("Calculation is carried out using {0} random numbers".format(10 * (indx+1)))
print()

### using uniform distribution
lowb = 0; highb = 1.0
sumv = 0.0; sumsv = 0.0
for ii in range(numsample):
    y = np.random.uniform(lowb, highb);
    sumv += fx(y)
    sumsv += (fx(y))**2 
integral_val = (highb - lowb)*sumv/numsample
print("Value of the definite integral assuming a uniform distribution: {0:2.4f}".format(integral_val))
error_sigma = (highb - lowb) * np.sqrt((sumsv/numsample-(sumv/numsample)**2)) / np.sqrt(numsample)
print("Estimated standard deviation: {0:2.6f}".format(error_sigma))
print("Within 3 std: {0:2.3f}-{1:2.3f}". format(integral_val-3.0*error_sigma, integral_val+3.0*error_sigma))
print("Calculation is carried out using {0} random numbers".format(numsample))
print()

### analytical solution
I = lambda x: 4*np.arctan(x)
integral_val = I(1) - I(0)
print("Value of the definite integral from analytical solution: {0:2.4f}".format(integral_val))