Create a python script to run a month of MODIS-CALIOP collocated data

					
#!/usr/bin/env python

import numpy as np
import calendar
import os
import glob
import re

command = ("gfortran ModisCaliopCollocatedPixels.f90"
           " -I/hdf-4.2.4-snap4/hdf4/include -L/hdf-4.2.4-snap4/hdf4/lib"
           " -lmfhdf -ldf -ljpeg -lz"
           " -o ModisCaliopCollocatedPixels")

os.system(command)

year = 2008
month = 1

path_to_myd03 = '/Volumes/Data/Myd03/'
path_to_01km_clay = '/Volumes/Data/Caliop/01km_CLay/'
path_to_results = '/Users/marchant/Desktop/CollocatedModisCaliopPixels_db/'

caliop_re = re.compile("T(.*)-(.*)-(.*).hdf")

for day in np.arange(calendar.monthrange(int(year),int(month))[1]):
#for i in range(8,calendar.monthrange(int(year),int(month))[1],1):

	d = str(year) + '/' + str(year) + '_' + "%02d" % month + '_' + "%02d" % (day+1) + '/'

	path_to_myd03_day = path_to_myd03 + d
	path_to_01km_clay_day = path_to_01km_clay + d
	
	list_of_myd03_files = glob.glob(path_to_myd03_day + '/*hdf' )
	list_of_01km_clay_files = glob.glob(path_to_01km_clay_day + '/*hdf' )

	for modis_granule_name in list_of_myd03_files:

			#offset = len(path_to_myd03_day)
			#modis_granule_hour = int( modis_granule_name[offset+15:offset+17] )
			#modis_granule_minute = int( modis_granule_name[offset+17:offset+19] )

			hm = modis_granule_name.split('.')[2]

			modis_granule_hour = int(hm[0:2])
			modis_granule_minute = int(hm[2:4])
						
			#-----find caliop files -----#
			minimum = 1000
			find_one = 0
			for caliop_granule_name in list_of_01km_clay_files:
					
				#offset = len(path_to_01km_clay_day)
				#caliop_granule_hour = int(caliop_granule_name[offset+47:offset+49])
				#caliop_granule_minute = int(caliop_granule_name[offset+50:offset+52])
				#print caliop_granule_hour, caliop_granule_minute
			
				obj = caliop_re.search(caliop_granule_name)
			
				caliop_granule_hour = int(obj.group(1))
				caliop_granule_minute = int(obj.group(2))
							
				difference = modis_granule_hour * 60 + modis_granule_minute - \
				            (caliop_granule_hour * 60 + caliop_granule_minute)
		
				if (difference < minimum ) and ( difference > 0 ):
						minimum = difference
						tmp_02 = caliop_granule_name
						find_one = 1

				print modis_granule_name
				print modis_granule_hour, modis_granule_minute
				print caliop_granule_name
				print caliop_granule_hour, caliop_granule_minute

			#----- create collocated caliop-modis pixel file -----#
			
			if ( find_one == 1 ):
					command = 'ModisCaliopCollocatedPixels' + ' ' + \
					           tmp_02 + ' ' + modis_granule_name
					os.system(command)
					command = 'mv calmod_index_' + "%02d" % modis_granule_hour + \
					          "%02d" % modis_granule_minute + ' ' + path_to_results + d
					os.system(command)

			print 'MODIS granule done !'
			print '--------------------'
	print 'Day ' + str(day+1) + ' done'

Plot MODIS C6 MYD04 Cloud Pixel Distance Land Ocean with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar
import pprint

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD04_L2.A2015190.0010.006.2015190191844.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Cloud_Pixel_Distance_Land_Ocean')

attributes = data_selected_id.attributes()

data = data_selected_id.get()

#print data
#print data_selected_id.info()
#pprint.pprint( attributes )
print data.max()
#print data.shape

_FillValue = attributes['_FillValue']

#----------------------------------------------------------------------------------------#
# Plot data

data[ data == _FillValue ] = 0.0

data = data * attributes['scale_factor']

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0,0.0001,10,20,30,40,50,60,70,80,90,100]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')                  
                         
plt.xticks([0,500,1000,1500,2000,2500], 
           ['0','500','1000','1500','2000','2500'], rotation=0, fontsize=8 )
plt.yticks([1000,2000,3000,4000], 
           ['1000','2000','3000','4000'], 
           rotation=0, fontsize=8 )

title = 'Cloud Pixel Distance Land Ocean \n MYD04 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

cbar_bounds = [0.0,10,20,30,40,50,60,70,80,90,100]
cbar_ticks = [0.0,10,20,30,40,50,60,70,80,90,100] 

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=8)

plt.savefig("myd04_c6_Cloud_Pixel_Distance_Land_Ocean.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Plot MODIS C6 MYD04 Optical Depth Land And Ocean with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar
import pprint

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD04_L2.A2015190.0010.006.2015190191844.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Optical_Depth_Land_And_Ocean')

attributes = data_selected_id.attributes()

data = data_selected_id.get()

#print data
#print data_selected_id.info()
#pprint.pprint( attributes )
#print data.max()
#print data.shape

_FillValue = attributes['_FillValue']

#----------------------------------------------------------------------------------------#
# Plot data

data[ data == _FillValue ] = 0.0

data = data * attributes['scale_factor']

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0,0.0001,0.1,0.2,0.3,0.4,0.5]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')                  
                         
plt.xticks([0,25,50,75,100,135], 
           ['0','25','50','75','100','135'], rotation=0, fontsize=8 )
plt.yticks([0,50,100,150,203], 
           ['0','50','100','150','203'], 
           rotation=0, fontsize=8 )

title = 'Optical Depth Land And Ocean \n MYD04 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

cbar_bounds = [0.0,0.1,0.2,0.3,0.4,0.5]
cbar_ticks = [0.0,0.1,0.2,0.3,0.4,0.5] 

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=8)

plt.savefig("myd04_c6_Optical_Depth_Land_And_Ocean.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Download a MODIS MYD04 HDF granule file from ladsftp using python

					
#!/usr/bin/env python

from ftplib import FTP

import numpy as np
import os
import calendar
import pprint
import sys

InpArgLIst = sys.argv

#----------------------------------------------------------------------------------------#
# user inputs (case 1): date provided in the  script

year = 2015
month = 07
day = 9
hour = 00
minute = 10 

#----------------------------------------------------------------------------------------#
# user inputs (case 2): directly in command line

if len(InpArgLIst) == 6 :
	year = int( InpArgLIst[1] )
	month = int( InpArgLIst[2] )
	day = int( InpArgLIst[3] )
	hour = int( InpArgLIst[4] ) 
	minute = int( InpArgLIst[5] ) 

#----------------------------------------------------------------------------------------#
# Opening FTP connection

print( 'Opening FTP connection' )

ftp = FTP('ladsftp.nascom.nasa.gov')
#ftp = FTP('ftp://ladsweb.nascom.nasa.gov')
ftp.login('anonymous','')

day_of_year = 0 
for i in range(month-1):
    day_of_year = day_of_year + calendar.monthrange(year,i+1)[1]

day_of_year = day_of_year + day

day_of_year = "%03d" % day_of_year

ind = str(year) + day_of_year + '.' +  "%02d" % hour + "%02d" % minute

#----------------------------------------------------------------------------------------#
# get MYD04 C6 File

directory = 'allData/6/MYD04_L2/'+str(year)+'/'+day_of_year
ftp.cwd(directory)
print directory

yyy = []
ftp.retrlines('NLST', yyy.append)
#pprint.pprint(yyy)

matching = [s for s in yyy if ind in s]
print matching
print matching[0]

e = 1
n = 0 
while e == 1 and n < 10 :
	try:
		ftp.retrbinary('RETR ' + matching[0], open(matching[0], 'wb').write)
		e = 0 
	except:
		n = n + 1
		print("Error while downloading")

#----------------------------------------------------------------------------------------#
# Closing FTP connection

print 'Closing FTP connection'
ftp.close()

Plot MODIS C6 L2 Cloud Effective Radius (PCL and nonPCL) with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar

#----------------------------------------------------------------------------------------#
# Inputs

file_name = './inputs/myd06_granules_viz/MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

cer_id = file.select('Cloud_Effective_Radius')
cer_pcl_id = file.select('Cloud_Effective_Radius_PCL')

cer_data = cer_id.get()
cer_pcl_data = cer_pcl_id.get()

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Plot data

attributes_cer_id = cer_id.attributes()
attributes_cer_pcl_id = cer_pcl_id.attributes()

data_01 = cer_data * attributes_cer_id['scale_factor']
data_02 = cer_pcl_data * attributes_cer_pcl_id['scale_factor']

shape_data_01 = data_01.shape
for i in range(shape_data_01[0]):
	for j in range(shape_data_01[1]):
		if data_02[i,j] > 0: data_01[i,j] = data_02[i,j]

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-100,0,5,10,15,20,25,30,35,40,45,50,55,60]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data_01), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')

cbar_bounds = [0,5,10,15,20,25,30,35,40,45,50,55,60]
cbar_ticks = [0,5,10,15,20,25,30,35,40,45,50,55,60]     
               
cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=10)

title = 'Cloud Effective Radius (PCL + nonPCL) \n MYD06 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

plt.savefig("myd06_c6_cloud_effective_radius_PCL_and_nonPCL_baseline.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Plot MODIS C6 L2 Cloud Top Temperature 5km with python 3

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Cloud_Top_Temperature')

data = data_selected_id.get()

#----------------------------------------------------------------------------------------#
# Plot data

data = (data + 15000) * 0.01

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0,200,205,210,215,220,225,230,235,
          240,245,250,255,260,265,270,275,280,285,290,300]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')

cbar_bounds = [200,205,210,215,220,225,230,235,
               240,245,250,255,260,265,270,275,280,285,290,300]
cbar_ticks = [200,205,210,215,220,225,230,235,
              240,245,250,255,260,265,270,275,280,285,290,300]               
               
cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=10)

title = 'Cloud Top Temperature 5km \n MYD06 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

plt.savefig("myd06_c6_cloud_top_temperature_5km_baseline.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Gradient descent in fortran 90 (3D case)

					
program gradient_descent_3d

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i

real(kind=16) :: f,f_prime
real(pr) :: fp,h
real(pr) :: x1_0,x2_0,x3_0,z0,tmp_z0
real(pr) :: tmp_x1_0,tmp_x2_0,tmp_x3_0
integer :: nb_iter,nb_max_iter
real(pr) :: alpha,eps,cond 

alpha = 0.1 ! learning rate
nb_max_iter = 100 ! Nb max d'iteration
eps = 0.0001 ! stop condition

x1_0 = 0.0 ! start point
x2_0 = 0.5 
x3_0 = 0.0 
z0 = f(x1_0,x2_0,x3_0)

cond = eps + 10.0 ! start with cond greater than eps (assumption)
nb_iter = 0 
tmp_z0 = z0
do while( cond > eps .and. nb_iter < nb_max_iter)
	tmp_x1_0 = x1_0 - alpha * f_prime(1,x1_0,x2_0,x3_0)
	tmp_x2_0 = x2_0 - alpha * f_prime(2,x1_0,x2_0,x3_0)
	tmp_x3_0 = x3_0 - alpha * f_prime(3,x1_0,x2_0,x3_0)
	x1_0 = tmp_x1_0
	x2_0 = tmp_x2_0
	x3_0 = tmp_x3_0
	z0 = f(x1_0,x2_0,x3_0)
	nb_iter = nb_iter + 1
	cond = abs( tmp_z0 - z0 )
	tmp_z0 = z0	
	write(6,*) x1_0,x2_0,x3_0,cond
end do
	
end program gradient_descent_3d

real(kind=16) function f(x1,x2,x3)
implicit none
integer, parameter :: pr = selected_real_kind(15,3)
real(pr), intent(in) :: x1,x2,x3
real(pr) :: x1_sphere,x2_sphere,x3_sphere
real(pr) :: r1,r2,r3
x1_sphere = 1.5
x2_sphere = 2.5
x3_sphere = 0.5
r1 = x1 - x1_sphere
r2 = x2 - x2_sphere
r3 = x3 - x3_sphere
f = r1**2 + r2**2 + r3**2
end function f

real(kind=16) function f_prime(dpi,x1,x2,x3)
implicit none
integer, parameter :: pr = selected_real_kind(15,3)
integer, intent(in) :: dpi
real(pr), intent(in) :: x1,x2,x3
real(pr) :: x1_p,x2_p,x3_p
real(kind=16) :: f
real(pr) :: h
h = 0.001
x1_p = x1
x2_p = x2
x3_p = x3
if( dpi == 1 ) x1_p = x1_p + h
if( dpi == 2 ) x2_p = x2_p + h
if( dpi == 3 ) x3_p = x3_p + h
f_prime = ( f(x1_p,x2_p,x3_p) - f(x1,x2,x3) ) / h
end function f_prime

Gradient descent in fortran 90 (1D case)

					
program gradient_descent_1d

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i

real(kind=16) :: f,f_prime
real(pr) :: fp,h
real(pr) :: x0,y0,tmp_y 
integer :: nb_iter,nb_max_iter
real(pr) :: alpha,eps,cond 

alpha = 0.1 ! learning rate
nb_max_iter = 100 ! Nb max d'iteration
eps = 0.0001 ! stop condition

x0 = 1.5 ! start point
y0 = f(x0) 

cond = eps + 10.0 ! start with cond greater than eps (assumption)
nb_iter = 0 
tmp_y = y0
do while( cond > eps .and. nb_iter < nb_max_iter)
	x0 = x0 - alpha * f_prime(x0)
	y0 = f(x0)
	nb_iter = nb_iter + 1
	cond = abs( tmp_y - y0 )
	tmp_y = y0
	write(6,*) x0,y0,cond
end do
	
end program gradient_descent_1d

real(kind=16) function f(x)
implicit none
integer, parameter :: pr = selected_real_kind(15,3)
real(pr), intent(in) :: x
f = 3.0*x**2.0+2.0*x+1.0
end function f

real(kind=16) function f_prime(x)
implicit none
integer, parameter :: pr = selected_real_kind(15,3)
real(pr), intent(in) :: x
real(kind=16) :: f
real(pr) :: h
h = 0.001
f_prime = ( f(x+h) - f(x) ) / h
end function f_prime

Gaussian Processes for regression and classification (2d example) in python 2.7

					
# Author: Benjamin Marchant

from numpy.linalg import inv
from scipy import misc

import numpy as np
import matplotlib.pyplot as plt
import math

#----------------------------------------------------------------------------------------#
# Training data

x1_data = np.array([-8.0, -6.00, -7.0, -4.0, 2.0, 5.0, 1.0, 3.0, 7.0])
x2_data = np.array([-8.0, -3.0, 2.0, 4.0, 3.0, 7.0, -1.0, -4.0, -7.0])

y_data = np.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0])

sigma_n = 0.1

x1_min = min(x1_data) - 2.0
x1_max = max(x1_data) + 2.0

x2_min = min(x2_data) - 2.0
x2_max = max(x2_data) + 2.0

matplotlib_marker = []
for i in y_data:
	if i == -1.0: matplotlib_marker.append('o')
	if i == 1.0: matplotlib_marker.append('x')

plt.xlabel('x1')
plt.ylabel('x2')
plt.xlim(x1_min,x1_max)
plt.ylim(x2_min,x2_max)
plt.title('Training Data')

plt.grid(True)

for i in range(len(x1_data)):
	plt.scatter(x1_data[i], x2_data[i], marker=matplotlib_marker[i])

plt.savefig('gaussian_processes_2d_training_data.png', bbox_inches='tight')
#plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Find Gaussian processes (3) hyperparameters using gradient descent

def log_likelihood_function(gph_l1, gph_l2, gph_sigma_f):
	dim_x_data = len(x1_data)
	k = np.zeros((dim_x_data,dim_x_data))
	for i in range(dim_x_data):
		for j in range(i+1):
			dx1 = x1_data[i] - x1_data[j]
			dx2 = x2_data[i] - x2_data[j]
			coef = - 0.5 * ( dx1**2 / gph_l1**2 + dx2**2 / gph_l2**2 ) 
			k[i,j] = gph_sigma_f * gph_sigma_f * math.exp( coef )
			k[j,i] = k[i,j]
			if i == j:
				k[i,j] = k[i,j] + sigma_n * sigma_n 
	k_inv = inv(k)
	m1 = np.dot(k_inv,y_data)
	part1 = np.dot(y_data.T,m1)
	part2 = math.log(np.linalg.det(k))
	return -0.5 * part1 - 0.5 * part2 - dim_x_data / 2.0 * math.log(2*math.pi)

def partial_derivative(func, var=0, point=[]):
    args = point[:]
    def wraps(x):
        args[var] = x
        return func(*args)
    return misc.derivative(wraps, point[var], dx = 1e-6)

alpha = 0.1 # learning rate
nb_max_iter = 100 # Nb max of iterations
eps = 0.0001 # stop condition

gph_l1 = 2.5 # start point
gph_l2 = 2.5 # start point
gph_sigma_f = 3.0 # start point

ll_value = log_likelihood_function(gph_l1, gph_l2, gph_sigma_f)

cond = eps + 10.0 # start with cond greater than eps (assumption)
nb_iter = 0 
tmp_ll_value = ll_value
while cond > eps and nb_iter < nb_max_iter:
	
	tmp_gph_l1 = gph_l1 + \
	alpha * partial_derivative(log_likelihood_function, 0, [gph_l1, gph_l2, gph_sigma_f])
	
	tmp_gph_l2 = gph_l2 + \
	alpha * partial_derivative(log_likelihood_function, 1, [gph_l1, gph_l2, gph_sigma_f])
	
	tmp_gph_sigma_f = gph_sigma_f + \
	alpha * partial_derivative(log_likelihood_function, 2, [gph_l1, gph_l2, gph_sigma_f])
	
	gph_l1 = tmp_gph_l1
	gph_l2 = tmp_gph_l2
	gph_sigma_f = tmp_gph_sigma_f
	
	ll_value = log_likelihood_function(gph_l1, gph_l2, gph_sigma_f)
	
	nb_iter = nb_iter + 1
	
	cond = abs( tmp_ll_value - ll_value )
	
	tmp_ll_value = ll_value
	
	print 'gph_l1, gph_l2, gph_sigma_f, nb_iter, cond', \
	gph_l1, gph_l2, gph_sigma_f, nb_iter, cond

print "Hyperparameters found using gradient descent: "
print 'gph_l1, gph_l2, gph_sigma_f', gph_l1, gph_l2, gph_sigma_f

#----------------------------------------------------------------------------------------#
# Prediction 1 point  

dim_x1_data = len(x1_data)
k = np.zeros((dim_x1_data,dim_x1_data))
for i in range(dim_x1_data):
	for j in range(i+1):
		dx1 = x1_data[i] - x1_data[j]
		dx2 = x2_data[i] - x2_data[j]
		coef = - 0.5 * ( dx1**2 / gph_l1**2 + dx2**2 / gph_l2**2 ) 	
		k[i,j] = gph_sigma_f * gph_sigma_f * math.exp( coef )
		k[j,i] = k[i,j]
		if i == j:
			k[i,j] = k[i,j] + sigma_n * sigma_n 

x1_new = -7.0
x2_new = -5.0
k_new = np.zeros((dim_x1_data))
for i in range(dim_x1_data):
	dx1 = x1_new - x1_data[i]
	dx2 = x2_new - x2_data[i]	
	coef = - 0.5 * ( dx1**2 / gph_l1**2 + dx2**2 / gph_l2**2 ) 
	k_new[i] = gph_sigma_f * gph_sigma_f * math.exp( coef )

k_inv = inv(k)
m1 = np.dot(k_new,k_inv)
y_new = np.dot(m1,y_data)

print "Prediction for 1 point: "
print 'x1_new, x2_new, y_new', x1_new, x2_new, y_new

var_y = k[0,0] - k_new.dot(k_inv.dot(np.transpose(k_new)))
#print "var_y", var_y

x1 = np.append(x1_data,x1_new)
x2 = np.append(x2_data,x2_new)
y = np.append(y_data,y_new)

matplotlib_marker = []
for i in y:
	if i == -1.0: matplotlib_marker.append('o')
	if i == 1.0: matplotlib_marker.append('x')
	if i != -1.0 and i != 1.0: matplotlib_marker.append('d')

plt.xlim(x1_min, x1_max)
plt.ylim(x2_min, x2_max)

plt.title('Prediction for a new point')

plt.xlabel('x1')
plt.ylabel('x2')

plt.grid(True)

print len(x1), len(x2), len(matplotlib_marker) 

for i in range(len(x1)):
	plt.scatter(x1[i], x2[i], marker=matplotlib_marker[i])

plt.savefig('gaussian_processes_2d_prediction_1_point.png', bbox_inches='tight')
#plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Inference regular gridded data

x1 = np.arange(-10.0,10.0,0.1)
x2 = np.arange(-10.0,10.0,0.1)

dim_x1 = x1.shape[0]
dim_x2 = x2.shape[0]
dim_x = dim_x1 * dim_x2

Z = np.zeros((dim_x))
Z_var = np.zeros((dim_x))

i1_cp = 0
i2_cp = 0
for i in range(dim_x):
	x1_new = x1[i1_cp]
	x2_new = x2[i2_cp]
	k_new = np.zeros((dim_x1_data))
	for j in range(dim_x1_data):
		dx1 = x1_new - x1_data[j] 
		dx2 = x2_new - x2_data[j]
		coef = - 0.5 * ( dx1**2 / gph_l1**2 + dx2**2 / gph_l2**2 ) 
		k_new[j] = gph_sigma_f * gph_sigma_f * math.exp( coef )
	k_inv = inv(k)
	m1 = np.dot(k_new,k_inv)
	y_new = np.dot(m1,y_data)
	Z[i] = y_new
	Z_var[i] = k[0,0] - k_new.dot(k_inv.dot(np.transpose(k_new)))
	i2_cp = i2_cp + 1
	if i2_cp == dim_x2:
		i1_cp = i1_cp + 1
		i2_cp = 0
	if i1_cp == dim_x1:
		i1_cp = 0

Z = np.reshape(Z, (dim_x1,dim_x2)) 
Z_var = np.reshape(Z_var, (dim_x1,dim_x2)) 

plt.imshow(Z.T, interpolation='bilinear', origin='lower')	
plt.colorbar()
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Gaussian processes 2d regression')
plt.grid(True)
plt.savefig('gaussian_process_2d_gridded_data.png')
plt.close()

plt.imshow(Z_var.T, interpolation='bilinear', origin='lower')	
plt.colorbar()
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Gaussian processes 2d regression (variance)')
plt.grid(True)
plt.savefig('gaussian_process_2d_gridded_data_var.png')
plt.close()

#----------------------------------------------------------------------------------------#
# Classification

for x in np.nditer(Z, op_flags=['readwrite']):
	x[...] = 1.0 / ( 1.0 + math.exp(-x))

plt.imshow(Z.T, interpolation='bilinear', origin='lower',vmin=0,vmax=1)	
plt.colorbar()
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Gaussian processes 2d for classification')
plt.grid(True)
plt.savefig('gaussian_process_2d_gridded_data_normalized.png')
plt.close()

CS = plt.contour(x1, x2, Z.T, origin='lower')
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Gaussian processes 2d for classification (contour plot)')
plt.grid(True)
plt.savefig('gaussian_process_2d_gridded_data_normalized_contour_plot.png')
plt.close()

Scattering geometry illustration for a non-spherical particle with metapost

					
verbatimtex
%&latex
\documentclass{article}
\usepackage[english,frenchb]{babel}
\begin{document}
etex

outputformat := "png";
outputtemplate := "%j-%c.%o";

input 3danim
output_res:=72;
drawing_scale:=1cm;

set_plane_(projection_plane)(point_null,vec_I,vec_K);
new_plane_(horizontal_plane)(point_null,vec_I,vec_J);

numeric u,r,dl,const,tetarot;
u=2;
r=0.25u;
dl=0.5u;
const= 0.866025403;
tetarot=-45;

vardef def_C(expr inst)=
  new_obj_points(inst,27); 
  set_C_points(inst);
enddef;

vardef set_C_points(expr inst)=

set_point(1)(cosd(tetarot)*(r)-sind(tetarot)*(dl),0,sind(tetarot)*(r)+cosd(tetarot)*(dl));
set_point(2)(cosd(tetarot)*(0.5r)-sind(tetarot)*(dl),r*const,sind(tetarot)*(0.5r)+cosd(tetarot)*(dl));
set_point(3)(cosd(tetarot)*(-0.5r)-sind(tetarot)*(dl),r*const,sind(tetarot)*(-0.5r)+cosd(tetarot)*(dl));
set_point(4)(cosd(tetarot)*(-r)-sind(tetarot)*(dl),0,sind(tetarot)*(-r)+cosd(tetarot)*(dl));
set_point(5)(cosd(tetarot)*(-0.5r)-sind(tetarot)*(dl),-r*const,sind(tetarot)*(-0.5r)+cosd(tetarot)*(dl));
set_point(6)(cosd(tetarot)*(0.5r)-sind(tetarot)*(dl),-r*const,sind(tetarot)*(0.5r)+cosd(tetarot)*(dl));

set_point(7)(cosd(tetarot)*(r)-sind(tetarot)*(-dl),0,sind(tetarot)*(r)+cosd(tetarot)*(-dl));
set_point(8)(cosd(tetarot)*(0.5r)-sind(tetarot)*(-dl),r*const,sind(tetarot)*(0.5r)+cosd(tetarot)*(-dl));
set_point(9)(cosd(tetarot)*(-0.5r)-sind(tetarot)*(-dl),r*const,sind(tetarot)*(-0.5r)+cosd(tetarot)*(-dl));
set_point(10)(cosd(tetarot)*(-r)-sind(tetarot)*(-dl),0,sind(tetarot)*(-r)+cosd(tetarot)*(-dl));
set_point(11)(cosd(tetarot)*(-0.5r)-sind(tetarot)*(-dl),-r*const,sind(tetarot)*(-0.5r)+cosd(tetarot)*(-dl));
set_point(12)(cosd(tetarot)*(0.5r)-sind(tetarot)*(-dl),-r*const,sind(tetarot)*(0.5r)+cosd(tetarot)*(-dl));

set_point(13)(0u,0u,0u);
set_point(14)(-3u,0u,0u);
set_point(15)(4u,4u,0u);
set_point(16)(0u,0u,3u);

set_point(17)(cosd(tetarot)*(0)-sind(tetarot)*(3u),0u,sind(tetarot)*(0)+cosd(tetarot)*(3u));
set_point(18)(cosd(tetarot)*(0)-sind(tetarot)*(dl),0u,sind(tetarot)*(0)+cosd(tetarot)*(dl));

set_point(19)(cosd(tetarot)*(0)-sind(30)*(2.5u),0u,sind(tetarot)*(0)+cosd(30)*(2.5u));
set_point(20)(cosd(tetarot)*(0)-sind(130)*(1.9u),0u,sind(tetarot)*(0)+cosd(130)*(1.9u));
set_point(21)(cosd(tetarot)*(0)-sind(-80)*(2u),0u,sind(tetarot)*(0)+cosd(-80)*(2u));

set_point(22)(cosd(-60)*(0)-sind(-60)*(3u),0u,sind(-60)*(0)+cosd(-60)*(3u));

set_point(23)(0u,-0.1u,0.65u);
set_point(24)(0u,-0.8u,0.8u);
set_point(25)(0u,-2u,1.3u);
set_point(26)(0u,1.6u,2.6u);
set_point(27)(0u,-1.5u,1u);

enddef;

vardef draw_C(expr inst)=

draw_line(1,2);
draw_line(2,3);
draw_line(3,4);
draw_line(4,5);
draw_line(5,6);
draw_line(6,1);

draw_line(7,8);
draw_line(8,9);
draw_line(9,10);

drawoptions(dashed evenly);

draw_line(10,11);
draw_line(11,12);
draw_line(12,7);

draw_line(13,18);

drawoptions();

draw_line(1,7);
draw_line(2,8);
draw_line(3,9);
draw_line(4,10);

draw_line(13,14);
draw_line(13,15);
draw_line(13,16);

draw_line(18,17);

draw_arrow(19,13);
draw_line(13,20);
draw_line(13,21);

draw_arrow(13,22);

%drawarrow (0.25z22){dir(60)}..(0.20z19){dir(-40)};
drawarrow (0.35z17){dir(40)}..(0.40z19){dir(-40)};
drawarrow (0.65z16){dir(0)}..(0.75z19){dir(-40)};
drawarrow (0.65z16){dir(180)}..(0.75z22){dir(230)};
drawarrow (0.50z16){dir(180)}..(0.50z17){dir(230)};


label_obj.rt(btex $Y$ etex,14);
label_obj.llft(btex $X$ etex,15);
label_obj.top(btex \begin{tabular}{c p{4.7cm}} 
Zenith \\
Z
\end{tabular} etex,16);

label_obj.top(btex \begin{tabular}{c p{4.7cm}} 
Particle Orientation \\
$Z_s$
\end{tabular} etex,17);

label_obj.top(btex \begin{tabular}{c p{4.7cm}} 
Incident Beam \\
Z' 
\end{tabular}etex,19);

label_obj.ulft(btex Scattered Beam etex,22);

label_obj.rt(btex $Y'$ etex,20);
label_obj.lft(btex $X'$ etex,21);

%label_obj.lft(btex $\pi-\Theta$ etex,23);
label_obj.lft(btex $\alpha'$ etex,24);
label_obj.lft(btex $\theta'$ etex,25);
label_obj.lft(btex $\theta$ etex,26);
label_obj.lft(btex $\alpha$ etex,27);

enddef;

assign_obj("cube","C");
%point_of_view_obj("cube",4,Obs_phi);

oblique_projection(horizontal_plane)(pnt_obj("cube",13),CABINET,45,45);

beginfig(100);

Obs_dist:=2;
draw_obj("cube"); 

endfig;

end.