##  octopus epidermis pattern mimicry mechanisms model 2021.1.1
##  Takeshi Ishida
import pygame
import random
import math
from pygame.locals import *
import sys
import openpyxl as px

# -----------------------------------------------------
# Function that multiplies a filter to calculate a Turing pattern 
def rule(cells,filter1,filter2):

  # Generate a two-dimensional list that holds the following states 
  tcells = [[0] * (meshNum) for i in range(meshNum)]
  # Filter convolution operation 
  for y in range(0, meshNum):
    for x in range(0, meshNum):
      sum1 =0
      sum2 =0
      for i in range(0, meshFilter):
        for j in range(0, meshFilter):
           sum1 += cells[(y+i-(meshFilter//2)+meshNum)%meshNum][(x+j-(meshFilter//2)+meshNum)%meshNum] * filter1[i][j]
           sum2 += cells[(y+i-(meshFilter//2)+meshNum)%meshNum][(x+j-(meshFilter//2)+meshNum)%meshNum] * filter2[i][j]

      # Applying transition rules 
      if (sum2*(1.0-w)-sum1*w)>0: tcells[y][x] = 1 
      if (sum2*(1.0-w)-sum1*w)<0: tcells[y][x] = 0 
      if (sum2*(1.0-w)-sum1*w)==0: tcells[y][x] = cells[y][x] 
 
  return tcells # Returns a next-step 2D list 

# -----------------------------------------------------
# Function that inverses a parameter from a Turing pattern 
def recal(cells,rr):
  # Prepare directions to access adjacent cells
  dir = ((-1, -1), ( 0, -1), ( 1, -1),
         (-1,  0),           ( 1,  0),
         (-1,  1), ( 0,  1), ( 1,  1))
		 
  # Generate a two-dimensional list that holds the w value of the boundary 
  www = [[0] * (meshNum) for i in range(meshNum)]
  # Generate a 2D list that holds the boundary status 
  blackborder = [[0] * (meshNum) for i in range(meshNum)]
  whiteborder = [[0] * (meshNum) for i in range(meshNum)]

  # Generate a two-dimensional list of filters for inverse operations 
  filterA = [[0] * (meshFilter) for i in range(meshFilter)]  #Outer neighborhood filter 
  filterB = [[0] * (meshFilter) for i in range(meshFilter)]  #Inner neighborhood filter  
  filterlen2 = [[0] * (meshFilter) for i in range(meshFilter)]  # Distance from the center of the filter to each cell 

  #　Creating filters 
  for y in range(0, meshFilter):
    for x in range(0,meshFilter):
      i = (meshFilter // 2)*(-1)+y
      j = (meshFilter // 2)*(-1)+x
	
      filterlen2[y][x] = math.sqrt(i*i + j*j)

      if filterlen2[y][x] <= rr:
        filterA[y][x] = 1
      if filterlen2[y][x] <= rr/2:
        filterB[y][x] = 1
        
  # Calculation of w value at black and white boundary 
  sumB =0
  sumW =0
  countB =0
  countW =0

  for y in range(0, meshNum):
    for x in range(0, meshNum):
      # Boundary judgment 
      blackborder[y][x]=0 
      whiteborder[y][x]=0
      if cells[y][x]==1:
        c = 0
        for d in dir: # Judge the 1 (black) border. If there is even one 0 next to 1, a black border 
          if cells[((y + d[1])+meshNum)%meshNum][((x + d[0])+meshNum)%meshNum] == 1: c += 1
        if c<8 : blackborder[y][x] = 1

      if cells[y][x]==0:
        c = 0
        for d in dir: # # Judge the 0 (white) border. If there is even one 1 next to 0, a white border
          if cells[((y + d[1])+meshNum)%meshNum][((x + d[0])+meshNum)%meshNum] == 0: c += 1
        if c<8 : whiteborder[y][x] = 1

      # Inverse operation of w value (ww) 
      sum11 =0
      sum22 =0
      for i in range(0, meshFilter):
        for j in range(0, meshFilter):
           sum11 += cells[(y+i-(meshFilter//2)+meshNum)%meshNum][(x+j-(meshFilter//2)+meshNum)%meshNum] * filterA[i][j]
           sum22 += cells[(y+i-(meshFilter//2)+meshNum)%meshNum][(x+j-(meshFilter//2)+meshNum)%meshNum] * filterB[i][j]
          
      if sum11 != 0 and sum22 != 0:
        www[y][x] = sum22/(sum11+sum22)

        if blackborder[y][x]==1 : 
          sumB += www[y][x]
          countB += 1	
#        print(rr,y,x,www[y][x])
        if whiteborder[y][x]==1 : 
          sumW += www[y][x]
          countW += 1			
#        print(rr,y,x,www[y][x])
      else : www[y][x] =0
  ww =(sumB+sumW)/(countB+countW)	

  # Index value calculation 		
  countBu =0
  countWu =0
  ss =0

  for y in range(0, meshNum):
    for x in range(0, meshNum):
      if blackborder[y][x]==1 and www[y][x]>=ww : 
            countBu += 1
      if whiteborder[y][x]==1 and www[y][x]<=ww : 
            countWu += 1			

      if blackborder[y][x]==1 : 
            ss = ss + (www[y][x]-ww)**2			
      if whiteborder[y][x]==1 : 
            ss = ss + (www[y][x]-ww)**2			
	  
  # calculation of Index value 		
  aa = countBu/countB + countWu/countW  
  ss = math.sqrt(ss/(countBu+countWu))
  print(rr,aa,ss,ww)
    
  return ww,aa,ss # Returns an estimate values 
  
# -------------------------------
#    Main program 
# -------------------------------
  
# Program initialization 
pygame.init()

meshNum = 40       #Number of vertical and horizontal cells 
meshPich = 20      #Drawing interval of cell  
meshFilter = 17    #Set the filter size and an odd number less than meshNum 
r = 3.0            #The radius of the outer neighborhood. Set less than the size of the filter 
w = 0.25           #Shape parameter w to create an image to reproduce

shigma =0

# Generate drawing screen 
screen = pygame.display.set_mode((((meshNum)*meshPich)*2+meshPich+100, ((meshNum)*meshPich)+100))

# Generate button  
button  = pygame.Rect(20, 20, 70, 50)   # creates a rect object
button2 = pygame.Rect(100, 20, 70, 50)  # creates a rect object
button3 = pygame.Rect(180, 20, 70, 50)  # creates a rect object

font = pygame.font.Font(None, 24)
text1 = font.render("Make", True, (0,0,0))   # Button to generate the original image of Turing pattern 
text2 = font.render("Calc", True, (0,0,0))   # Button to back-calculate parameters from Turing pattern image 
text3 = font.render("Reverse", True, (0,0,0)) # Button to reproduce an image from the inverse calculation parameters 

pygame.display.set_caption("Octopus Model 2020")    # Create title 

# Generate a two-dimensional list that holds the state of the cells 
cells = [[0] * (meshNum) for i in range(meshNum)]
xcells = [[0] * (meshNum) for i in range(meshNum)]

# Set a random number (0,1) in the 2D list 
for y in range(0, meshNum):
  for x in range(0, meshNum):
    cells[y][x] = random.randrange(2)
    xcells[y][x] = random.randrange(2)

# Reading images from Excel files (when using a pattern without Turing pattern) 
#wb1 = px.load_workbook('C:/SeaBed/SeaBed40x40-5.xlsx', data_only=True)
#wb1 = px.load_workbook('C:/SeaBed/chekerbord10x10hensoku.xlsx', data_only=True)
#ws1 = wb1.active

#for i in range(0,meshNum):
#  for j in range(0,meshNum):
#      cells[i][j]=(ws1.cell(column=(j+1), row=(i+1)).value)


# Generate a two-dimensional list of filter layers 
filter1 = [[0] * (meshFilter) for i in range(meshFilter)]  #Outer neighborhood filter
filter2 = [[0] * (meshFilter) for i in range(meshFilter)]  #Inner neighborhood filter 
filterlen = [[0] * (meshFilter) for i in range(meshFilter)]  # Distance from the center of the filter to each cell 

# Generate a two-dimensional list of filter layers for inverse operations 
filterA = [[0] * (meshFilter) for i in range(meshFilter)]  #Outer neighborhood filter
filterB = [[0] * (meshFilter) for i in range(meshFilter)]  #Inner neighborhood filter   


#　Creating filterｓ 
for y in range(0, meshFilter):
  for x in range(0,meshFilter):
    i = (meshFilter // 2)*(-1)+y
    j = (meshFilter // 2)*(-1)+x
	
    filterlen[y][x] = math.sqrt(i*i + j*j)

    if filterlen[y][x] <= r:
      filter1[y][x] = 1
    if filterlen[y][x] <= r/2:
      filter2[y][x] = 1
	
#print(filter1)	
#print(filter2)	
	
run = 0 # Execution flag 
gen = 1 # 

while True: # Infinite loop to run the simulation 
  screen.fill((255, 255, 255)) # Make the background white 

  # display of buttons   
  pygame.draw.rect(screen, (255, 0, 0), button)
  pygame.draw.rect(screen, (0, 255, 0), button2)
  pygame.draw.rect(screen, (0, 0, 255), button3)

  screen.blit(text1, (25, 35))
  screen.blit(text2, (105,35))
  screen.blit(text3, (185,35))

  if run == 1: 
    gen += 1

    cells = rule(cells,filter1,filter2)    #Turing pattern calculation 

    pygame.time.wait(150) # Stop the program for 150 milliseconds 


  # Draw cell states 
  for y in range(0, meshNum):
    for x in range(0, meshNum):
      if cells[y][x] == 1:
        rect = Rect(x * meshPich+80, y * meshPich+80, meshPich, meshPich)
        pygame.draw.rect(screen, (128, 128, 128), rect)
      if xcells[y][x] == 1:
        rect = Rect(x * meshPich+meshNum*meshPich+meshPich+80, y * meshPich+80, meshPich, meshPich)
        pygame.draw.rect(screen, (128, 128, 128), rect)
		
  pygame.draw.rect(screen, (0, 0, 0), ( 80, 80, meshNum * meshPich, meshNum * meshPich ),2)
  pygame.draw.rect(screen, (0, 0, 0), ( meshNum * meshPich+meshPich+80, 80, meshNum * meshPich, meshNum * meshPich ),2)
  
  # Process mouse input 
  for e in pygame.event.get():
     if e.type == pygame.MOUSEBUTTONDOWN:

        if button.collidepoint(e.pos):
          # Create a Turing patterns, Show pattern on the left side of the screen 
          run = 1 - run # Flip the run flag 

        if button2.collidepoint(e.pos):
          # Inverse calculation of Turing pattern parameters 
          xamax =0
          xrmax =0
          xwmax =0
          xsmax =10000
    
          for i in range(1, 11):
             xr = 8.0-i*0.5 
             xw1,xa1,shigma = recal(cells,xr)
             if xamax<xa1 :
#             if xsmax>shigma :
               xamax = xa1
               xrmax = xr
               xwmax = xw1
               xsmax = shigma
          print(xrmax,xwmax,xamax,xsmax)
          
        if button3.collidepoint(e.pos):
          # Reproduce the pattern on the right part of the screen from the inverse calculated value 
          #run=1
          for y in range(0, meshFilter):
            for x in range(0,meshFilter):
               if filterlen[y][x] <= xrmax:
                  filter1[y][x] = 1
               if filterlen[y][x] <= xrmax/2:
                  filter2[y][x] = 1
          w=xwmax
          for y in range(0, 10):
             xcells = rule(xcells,filter1,filter2)    #Turing pattern calculation 
          
     if e.type == QUIT:
        pygame.quit()
        sys.exit()
      
  s1 = "running" if run == 1 else "setting" 
  s2 = "generation : " + str(gen)
  text = font.render(s1 + "   " + s2, True, (0, 0, 0))
  screen.blit(text, (1, 1)) 
  pygame.display.update()