Contact analysis. about gammaray HOT 2 CLOSED

´import sgems 
import numpy as np 
import scipy 
import matplotlib.pyplot as plt
import pylab as py


filename = "C:/Database/Recursos/JAMARI/3D/Misc/SGeMS/Contact analysis/results.dat" #INFORME O CAMINHO DO ARQUIVO DE SAÍDA
fid = file(filename,"w")	#command to write 

distlag = 0.5 #ENTRE COM O VALOR DO PRIMEIRO LAG. Se o lag for 10, o primeiro lag calcula as distância entre amostras de domínios diferentes que sejam menores ou iguais a 10
nlags = 15 #NÚMERO DE LAGS A SER CONSIDERADO. Os lags seguintes são múltiplos do primeiro lag

gridname = "SJA Dataset - LTBX" 
x = sgems.get_property(gridname,"_X_") #store coordinate X into array x 
y = sgems.get_property(gridname,"_Y_") #store coordinate Y into array y 
z = sgems.get_property(gridname,"_Z_") #store coordinate Z into array z 
dom = sgems.get_property(gridname,"LITN") #store domain into variable dom /ENTRE COM O NOME DO CARIMBO
var = sgems.get_property(gridname,"FEPRO") #store contents into variable var /ENTRE COM O NOME DO TEOR
a = np.vstack([x,y,z,var,dom]) 
array = a.transpose()

DomUm = 3 #INFORME O CÓDIGO DO PRIMEIRO DOMÍNIO
DomDois = 4 #INFORME O CÓDIGO DO SEGUNDO DOMÍNIO


dist= [] 
sum = []
sdois = []
for i in range(len(x)): 
	for j in range(i+1,len(x),1): 
		if dom[i] <> dom[j] and var[i] > -999 and var[j] > -999: 
			distancia = ((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2)**0.5	
			dist.append(distancia)
			sum.append(i)
			sdois.append(j)		
			

dist, sum,sdois = zip(*sorted(zip(dist,sum,sdois)))	 

lag = np.zeros(nlags)
minlag = np.zeros(nlags)
maxlag = np.zeros(nlags)
 
contum = np.zeros(nlags) 
contdois = np.zeros(nlags) 
contcont = np.zeros(nlags) 
MedUm = np.zeros(nlags)
MedDois = np.zeros(nlags)
DistMed = np.zeros(nlags)
DistRound = np.zeros(nlags)
AddDomUm = np.zeros(nlags) 
AddDomDois = np.zeros(nlags) 
SomaDist = np.zeros(nlags) 
MarkUm = np.zeros(len(var)) 
MarkDois = np.zeros(len(var))



for i in range(nlags): 
	lag[i] = ((i+1)*distlag)




for k in range (len(lag)): 
		
	

	for z in range(len(dist)):
		
			
		if dist[z] <= lag[k] and dist[z] <> -99 and MarkDois[sum[z]]==0 and MarkDois[sdois[z]]==0: #if distance equal or less the selected lag then... 
				
				if MarkUm[sum[z]] == 0 and MarkUm[sdois[z]] == 1: #primeiro caso, amostra i não está carimbada e amostra j está carimbada 
					MarkUm[sum[z]] = 1 
					contcont[k] = contcont[k] + 1 
					SomaDist[k] = SomaDist[k] + dist[z] 

					
#Amostra não carimbada está em qual domínio ? 1 ou 2 ? 
					if array[sum[z],4] == DomUm: # 'se amostra i está no domínio 1 então 
						
						AddDomUm[k] = AddDomUm[k] + array[sum[z], 3]# 'amostra i é adicionada à soma do domínio 1 
						AddDomDois[k] = AddDomDois[k]# ' soma do domínio dois permanece constante
						contum[k] = contum[k] + 1 #' contador de amostras do domínio 1 recebe mais uma unidade 
						contdois[k] = contdois[k] #' contador do domínio 2 permanece constante 
					else: # senão amostra j está no domínio dois 
						
						AddDomUm[k] = AddDomUm[k] #'soma do domínio 1 permanece constante 
						AddDomDois[k] = AddDomDois[k] + array[sum[z], 3] #' amostra i é adicionada à soma do domínio 2 
						contum[k] = contum[k] #' contador do domínio 1 permanece constante 
						contdois[k] = contdois[k] + 1 #' contador do dómínio 2 é acrescentado uma unidade 
				elif  MarkUm[sum[z]] == 1 and MarkUm[sdois[z]] == 0: 
# 'conta distancias calculadas 
					contcont[k] = contcont[k] + 1 
					SomaDist[k] = SomaDist[k] + dist[z] 
					MarkUm[sdois[z]] = 1 
#'Em qual domínio está amostra j? 
					if array[sdois[z], 4] == DomUm: # 'se ela está no domínio 1 
						contum[k] = contum[k] + 1 
						contdois[k] = contdois[k] 
						AddDomUm[k] = AddDomUm[k] + array[sdois[z], 3] 
						AddDomDois[k] = AddDomDois[k] 
					else:# 'senão ela está no domínio 2 
						contum[k] = contum[k] 
						contdois[k] = contdois[k] + 1 
						AddDomUm[k] = AddDomUm[k] 
						AddDomDois[k] = AddDomDois[k] + array[sdois[z], 3] 
				elif MarkUm[sum[z]]== 0 and MarkUm[sdois[z]] == 0: 
					SomaDist[k] = SomaDist[k] + dist[z]
					contcont[k] = contcont[k] + 1 
					contum[k] = contum[k] + 1 
					contdois [k]= contdois[k] + 1 
										
					MarkUm[sum[z]] = 1 
					MarkUm[sdois[z]] = 1 
					if array[sum[z],4] == DomUm: 
						AddDomUm[k] = AddDomUm[k] + array[sum[z],3] 
						AddDomDois[k] = AddDomDois[k] + array[sdois[z], 3] 
					else: 
						AddDomUm[k] = AddDomUm[k] + array[sdois[z], 3] 
						AddDomDois[k] = AddDomDois[k] + array[sum[z], 3]
				else:
					SomaDist[k] = SomaDist[k]
					contcont[k] = contcont[k]
					contum[k] = contum[k]
					contdois [k]= contdois[k]
	
	
	
	for i in range(len(MarkUm)):
		MarkDois[i] = MarkUm[i]			
			
				
	if contcont[k] == 0:
		contcont[k] = -99
		DistMed[k] = -99
		contcont[k] = -99
	else:		
		contcont[k] = contcont[k]
		DistMed[k] = SomaDist[k]/(2*contcont[k])
		DistRound[k] = round(DistMed[k],2)

	if contum[k] == 0:
		MedUm[k] = -99
	else:
		MedUm[k] = AddDomUm[k]/contum[k]

	if contdois[k] == 0:
		MedDois[k] = -99
	else:
		MedDois[k] = AddDomDois[k]/contdois[k]		
		

	

separator = "	"	#indica que o separador utilizado será espaço
for k in range(nlags):
	if contcont[k] <> -99: 
		fid.write(str(MedUm[k])+separator+str(MedDois[k])+ separator+str(DistRound[k]) + separator+str(contum[k])+separator+str(contdois[k]) + "\n")
		
fid.close()





HighestX = 0
HighestYUm = 0
HighestYDois = 0
HighestY = 0
LowestX = 0
LowestY = 0

for p in range(nlags):
	if MedUm[p] > HighestYUm:
		HighestYUm = MedUm[p]

for p in range(nlags):
	if MedDois[p] > HighestYDois:
		HighestYDois = MedDois[p]

if HighestYDois > HighestYUm:
	HighestY = HighestYDois
else:
	HighestY = HighestYUm

for k in range(nlags):
	if DistRound[k] > HighestX:
		HighestX = DistRound[k]



for k in range(nlags):
	if DistMed[k] <> -99 and MedUm[k]<>-99 and MedDois[k]<>-99:
		

		
		py.plot(-DistMed[k],MedUm[k],'k^')
		py.plot(DistMed[k],MedDois[k],'bo')


py.ylim([0,(HighestY*1.05)])
py.xlabel("Distancia do contato")
py.ylabel("FEPRO")
py.legend(("LT","BX"),loc = "lower right")




plt.show()

import sgems

# Input parameters
gridname = 'grid'  # name of the dataset
grade_name = 'trend +X-'  # name of the property of interest
domain_name = 'lito'  # name of the property that contains the domains
tmin = -98

# inform the lag and the number of lags
lag = 1.0
nlags = 10

# inform the integer of the first and second domain
Domain_One = 1
Domain_Two = 2

X = sgems.get_property(gridname, '_X_')
Y = sgems.get_property(gridname, '_Y_')
Z = sgems.get_property(gridname, '_Z_')
Domain = sgems.get_property(gridname, domain_name)
Grade = sgems.get_property(gridname, grade_name)

# Inform output file
output_file = 'Contact_analysis.txt'


print('Starting script contact analysis')



x_graph_min = 0
x_graph_max = -20

y_graph_min = 0
y_graph_max = -20

max_dist = []
min_dist = []
for i in range(nlags):
    max_dist.append((i + 1) * lag)
    min_dist.append(i * lag)


def select_samples_within_domain(column_to_be_extracted, column_domain, column_grade, domain, tmin):
    column_extracted = []
    EPSLON = 0.000001
    for i in range(len(column_to_be_extracted)):
        if (abs(column_domain[i] - domain) < EPSLON and column_grade[i] > tmin and column_domain[i] > tmin):
            column_extracted.append(column_to_be_extracted[i])

    return column_extracted


X_one = select_samples_within_domain(X, Domain, Grade, Domain_One, tmin)
Y_one = select_samples_within_domain(Y, Domain, Grade, Domain_One, tmin)
Z_one = select_samples_within_domain(Z, Domain, Grade, Domain_One, tmin)
Grade_one = select_samples_within_domain(Grade, Domain, Grade, Domain_One, tmin)

X_two = select_samples_within_domain(X, Domain, Grade, Domain_Two, tmin)
Y_two = select_samples_within_domain(Y, Domain, Grade, Domain_Two, tmin)
Z_two = select_samples_within_domain(Z, Domain, Grade, Domain_Two, tmin)
Grade_two = select_samples_within_domain(Grade, Domain, Grade, Domain_Two, tmin)

dist = []
index_one = []
index_two = []

for i in xrange(len(Grade_one)):

    for j in xrange(i + 1, len(Grade_two)):
        distance = ((X_one[i] - X_two[j]) ** 2 + (Y_one[i] - Y_two[j]) ** 2 + (Z_one[i] - Z_two[j]) ** 2) ** 0.5

        dist.append(distance)
        index_one.append(i)
        index_two.append(j)

matrix = zip(dist, index_one, index_two)
matrix.sort()

for i in range(len(matrix)):
    dist[i] = matrix[i][0]
    index_one[i] = matrix[i][1]
    index_two[i] = matrix[i][2]

# initialize variables
sum_one = []
count_one = []
avg_one = []
sum_two = []
count_two = []
avg_two = []
sum_dist = []
count_dist = []
avg_dist = []

for i in xrange(nlags):
    sum_one.append(0.00)
    count_one.append(0)
    avg_one.append(0.00)

    sum_two.append(0.00)
    count_two.append(0)
    avg_two.append(0.00)

    sum_dist.append(0.00)
    count_dist.append(0)
    avg_dist.append(0.00)

# initialize variable that flags whether the sample was used or not to calculate the mean
mark_one = []
for i in range(len(Grade_one)):
    mark_one.append(0)

mark_two = []
for i in range(len(Grade_two)):
    mark_two.append(0)

EPSLON = 0.0000001

for k in range(nlags):
    maximum = max_dist[k]
    minimum = min_dist[k]

    for i in range(len(dist)):

        if (dist[i] <= maximum and dist[i] > minimum):
            grade_one = Grade_one[index_one[i]]
            grade_two = Grade_two[index_two[i]]

            sum_dist[k] = sum_dist[k] + dist[i]
            count_dist[k] = count_dist[k] + 1

            # consider the sample of domain one for mean if it was not already flagged
            if (abs(mark_one[index_one[i]] - 0) < EPSLON):
                sum_one[k] = sum_one[k] + grade_one
                count_one[k] = count_one[k] + 1
                mark_one[index_one[i]] = 1

            if (abs(mark_two[index_two[i]] - 0) < EPSLON):
                
                sum_two[k] = sum_two[k] + grade_two
                count_two[k] = count_two[k] + 1
                mark_two[index_two[i]] = 1

for k in range(nlags):

    if (abs(count_one[k] - 0) < EPSLON):
        avg_one[k] = -99
    else:
        avg_one[k] = sum_one[k] / count_one[k]

    if (abs(count_two[k] - 0) < EPSLON):
        avg_two[k] = -99
    else:
        avg_two[k] = sum_two[k] / count_two[k]

    if (abs(count_dist[k] - 0) < EPSLON):
        avg_dist[k] = -99
    else:
        avg_dist[k] = sum_dist[k] / count_dist[k]

# Get X and Y points for graph domain one
x_points_one = []
y_points_one = []
for i in range(nlags):
    if (avg_one[i] > -98):
        x_points_one.append((-avg_dist[i]) / 2)
        y_points_one.append(avg_one[i])

# Get x and y points for graph domain two
x_points_two = []
y_points_two = []
for i in range(nlags):
    if (avg_two[i] > -98):
        x_points_two.append((avg_dist[i]) / 2)
        y_points_two.append(avg_two[i])

with open(output_file, 'w') as fout:
    fout.write(
        'Dist_from_contact_{}, Mean_Grade_{},  Dist_from_contact_{}, Mean_Grade_{}  \n'.format(Domain_One, Domain_One,
                                                                                               Domain_Two, Domain_Two))
    for i in range(nlags):
        fout.write('{:.4f}, {:.4f}, {:.4f}, {:.4f} \n'.format(-avg_dist[i] / 2.0, avg_one[i], avg_dist[i] / 2.0, avg_two[i]))


print('Finish script contact analysis')