# -*- coding:Utf-8 -*-
from numpy import *
from random import *
# Claude TOUZET 11 - 2007
# SOFM code (DNA chip ICA)

learning=1 # 1 for learning ; 0 for propagation
nb_iterations=10 # learning iterations : 100
nb_entree=18  # nb d'entrees sur chaque neurone
nb_neuron_som1D=150 # 1-D neurons of the self-organizing map
nb_neuron_som=nb_neuron_som1D*nb_neuron_som1D # 2-D neurons of the self-organizing map
nb_samples=5 # number of learning samples for the SOM : 20446
SOM = zeros(nb_neuron_som)    
try : 
    dnachip_reduite = open ('dnachip_reduite1.dat', 'r')
except :
    print "Fichier dnachip_reduite1.dat introuvable"
    
s = dnachip_reduite.read()
dnachip_reduite.close()

ss=s.split()
#print ss
donnees = ones((nb_samples, nb_entree+16))
for i in range (0, nb_samples) : 
    for j in range (0, nb_entree+16) :
        donnees[i][j] = float(ss[i*nb_entree + j])
        # print '(',i,',',j,')= ',donnees[i,j]

CX=donnees[:,0:18]
#HI=donnees[:,18:33]
BA=CX[0:nb_samples,:] # BA : base d'apprentissage
AA=[] # ordre aléatoire des exemples de la B

print BA

if (learning==1) : # initialisation aléatoire des poids de la SOFM
    WSOM=zeros((nb_neuron_som , nb_entree))
    for i in range (0, nb_neuron_som) :
        for j in range (0, nb_entree) :
            WSOM[i,j] = 5.0 + random()/10

print WSOM

# Self-organizing map  
mu1=0.4 
beta1=0.4 
if (learning==1) : #LEARNING
    for n in range (0, nb_iterations) : # iterations
        print 'iteration = ',n
        mu = mu1 - mu1 * (n/nb_iterations)
        beta = beta1 - beta1 * (n/nb_iterations)
        for b in range (0, nb_samples) :
            a= floor (random()*(nb_samples-1)-0.001) +1 
            AA=(AA, a) 
            for m in range (0, nb_neuron_som) : # compute distances 
                SOM[m]=0
                for k in range (0, nb_entree) :
                    SOM[m] = SOM[m] + abs(BA[a,k]-WSOM[m,k]) 
             
            index = SOM.argmin()  # winner
            m=index; # winner 
            for k in range (0, nb_entree) :
                    WSOM[m,k] = WSOM[m,k] + mu*(BA[a,k]-WSOM[m,k]) 
            
            # bord gauche
            if (((index-1)/nb_neuron_som1D) != ceil ((index-1)/nb_neuron_som1D)) :
                m=index-1 # neighbour gauche
                if (m>0) :
                    for k in range (0,nb_entree) :
                            WSOM[m,k] = WSOM[m,k] + beta*(BA[a,k]-WSOM[m,k]) 
             
            if (((index)/nb_neuron_som1D) != ceil ((index)/nb_neuron_som1D)) :
                m=index+1 # neighbour droit 
                if (m<(nb_neuron_som+1)) :
                    for k in range (0,nb_entree) :
                            WSOM[m,k] = WSOM[m,k] + beta*(BA[a,k]-WSOM[m,k])
                            
            if (((index-nb_neuron_som1D)/nb_neuron_som1D) != ceil ((index-nb_neuron_som1D)/nb_neuron_som1D)) :
                m=index-nb_neuron_som1D # neighbour haut 
                if (m>0) :
                    for k in range (0,nb_entree) :
                             WSOM[m,k] = WSOM[m,k] + beta*(BA[a,k]-WSOM[m,k])
 
            if (((index+nb_neuron_som1D)/nb_neuron_som1D) != ceil ((index+nb_neuron_som1D)/nb_neuron_som1D)) :
                m=index+nb_neuron_som1D # neighbour bas  
                if (m<(nb_neuron_som+1)) :
                    for k in range (0,nb_entree) :
                            WSOM[m,k] = WSOM[m,k] + beta*(BA[a,k]-WSOM[m,k])

''' 
if (learning==1) :
    # Store the WSOM
    genes=open('WSOM_final.dat', 'w')
    for i in range (0, nb_neuron_som) : 
        for j in range (0, nb_entree) :
            genes.write (WSOM(i,j))
            
    genes.close()
    
'''