# Note on genome representation - p_int is a vector of integers of length n_gene,encoding gene type. 
# p_float is an n_gene x 8 array encoding the real values for each gene (each row is a gene)
# Hence p_int[1] is the associated type of gene p_float[1,:]

function initialise_grn(p_float::Array{Float64,2},p_int::Vector{Int64},grn_config::FGRNConfig)

    n_rec = 0
    n_beh = 0
    n_reg = 0
    n_env = 0
    
    n = grn_config.protein_size

    for i in 1:size(p_int,1)

        if p_int[i] == 1
            n_rec += 1
        elseif p_int[i] == 2
            n_beh += 1
        elseif p_int[i] == 3
            n_reg += 1
        else
            n_env +=1
        end
    end

    
    grn = ( rec_coding = ones(Int,n,n,n_rec),
            beh_coding = ones(Int,n,n,n_beh),
            reg_coding = ones(Int,n,n,n_reg),
            env_coding = ones(Int,n,n,n_env),

            rec_prom = ones(Int,n,n,n_rec),
            beh_prom = ones(Int,n,n,n_beh),
            reg_prom = ones(Int,n,n,n_reg),
            env_prom = ones(Int,n,n,n_env),

            rec_cod_thresh = zeros(Float64,n_rec),
            rec_aff_thresh = zeros(Float64,n_rec),
            rec_id = ones(Int,n_rec),

            beh_cod_thresh  = zeros(Float64,n_beh),
            beh_aff_thresh = zeros(Float64,n_beh),
            beh_id = ones(Int,n_beh),

            reg_cod_thresh  = zeros(Float64,n_reg),
            reg_aff_thresh = zeros(Float64,n_reg),
            reg_id = ones(Int,n_reg),

            env_cod_thresh = zeros(Float64,n_env),
            env_aff_thresh = zeros(Float64,n_env),
            env_id = ones(Int,n_env),

            Ct = grn_config.Ct, # assume average affinity = 0.3 ish, average affinity th - 0.5 as uniform
            Cs = grn_config.Cs, # assume average affinity = 0.3 ish, average affinity th - 0.5 as uniform
            Cw = grn_config.Cw,
            Cp = grn_config.Cp,
            Ci = grn_config.Ci,
           )
    
    n_rec = 0
    n_beh = 0
    n_reg = 0
    n_env = 0

    for i in 1:size(p_int,1)

        coding = mandelbrot_protein(p_float[i,1:3],grn_config.mandel_thresh,grn_config.max_iter,n)
        promoter = mandelbrot_protein(p_float[i,4:6],grn_config.mandel_thresh,grn_config.max_iter,n)
        threshold = p_float[i,7]
        affinity = p_float[i,8]

        if p_int[i] == 1
            n_rec += 1
            grn.rec_coding[:,:,n_rec] = coding
            grn.rec_prom[:,:,n_rec] = promoter
            grn.rec_cod_thresh[n_rec] = threshold
            grn.rec_aff_thresh[n_rec] = affinity
            grn.rec_id[n_rec] = i

        elseif p_int[i] == 2
            n_beh += 1
            grn.beh_coding[:,:,n_beh] = coding
            grn.beh_prom[:,:,n_beh] = promoter
            grn.beh_cod_thresh[n_beh] = threshold
            grn.beh_aff_thresh[n_beh] = affinity
            grn.beh_id[n_beh] = i
            
        elseif p_int[i] == 3
            n_reg += 1
            grn.reg_coding[:,:,n_reg] = coding
            grn.reg_prom[:,:,n_reg] = promoter
            grn.reg_cod_thresh[n_reg] = threshold
            grn.reg_aff_thresh[n_reg] = affinity
            grn.reg_id[n_reg] = i
            
        else
            n_env += 1
            grn.env_coding[:,:,n_env] = coding
            grn.env_prom[:,:,n_env] = promoter
            grn.env_cod_thresh[n_env] = threshold
            grn.env_aff_thresh[n_env] = affinity
            grn.env_id[n_env] = i
        end
    end
    
    return grn
    
end

function calculate_affinity(protein_a::Array{Int,2}, protein_b::Array{Int,2})
    
    return convert(Float64,sum(@. abs(protein_a - protein_b))) / (size(protein_a,1)*size(protein_a,1)*255)
    
end

function bind_probability(affinity::Float64,affinity_th::Float64,config::FGRNConfig)
    
    return (1 + tanh((affinity_th - affinity - config.Ct)/config.Cs))/2
    
end

function update_rule(product_conc::Float64,original_conc::Float64,config::FGRNConfig)
    
    return product_conc*tanh((product_conc - 0.5)/config.Cw)/config.Ci - original_conc/config.Cp
    
end

function updateConcentrations!(dc::Array{Float64,1},grn::NamedTuple, concentrations::Array{Float64,1},grn_config::FGRNConfig)
    
    product = zeros(Int,grn_config.protein_size,grn_config.protein_size)
    product_conc = zeros(Float64,grn_config.protein_size,grn_config.protein_size)
    
    fill!(dc, zero(Float64))
    
    nProteins::Int = 0
    affinity = 0 
    
    @inbounds for j in 1:grn_config.protein_size
        @inbounds for i in 1:grn_config.protein_size
            
            for k in 1:size(grn.rec_coding,3)
                
                if (concentrations[grn.rec_id[k]] > 0.5*grn.rec_cod_thresh[k]) && (grn.rec_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.rec_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.rec_id[k]]
                    
                    nProteins+=1
                    
                end
            end

            for k in 1:size(grn.reg_coding,3)
                
                if (concentrations[grn.reg_id[k]] > 0.5*grn.reg_cod_thresh[k]) && (grn.reg_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.reg_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.reg_id[k]]
                    
                    nProteins+=1
                    
                end
            end
                                
            for k in 1:size(grn.env_coding,3)
                
                if (concentrations[grn.env_id[k]] > 0.5*grn.env_cod_thresh[k]) && (grn.env_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.env_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.env_id[k]]

                    nProteins+=1

                end

            end
        end
    end
    
    product_conc = sum(product_conc)/(grn_config.protein_size*grn_config.protein_size)
    
    if nProteins != 0
        
        @inbounds for k in 1:size(grn.rec_prom,3)
                
            dc[grn.rec_id[k]] = 0
            
        end
                            
        @inbounds for k in 1:size(grn.beh_prom,3)
            
            affinity = calculate_affinity(grn.beh_prom[:,:,k],product)
            
            if affinity > grn.beh_aff_thresh[k]
                
                dc[grn.beh_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.beh_id[k]],grn_config)
            
            else
                
                dc[grn.beh_id[k]] = 0  
                
            end
        end
        
        
        @inbounds for k in 1:size(grn.reg_prom,3)
            
            affinity = calculate_affinity(grn.reg_prom[:,:,k],product)
            
            if affinity > grn.reg_aff_thresh[k]
                
                dc[grn.reg_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.reg_id[k]],grn_config)
            
            else
                
                dc[grn.reg_id[k]] = 0
            end
        end
        
                                       
        @inbounds for k in 1:size(grn.env_prom,3)
            
            affinity = calculate_affinity(grn.env_prom[:,:,k],product)
            
            if affinity > grn.env_aff_thresh[k]
                
                dc[grn.env_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.env_id[k]],grn_config)
            
            else
                
                dc[grn.env_id[k]] = 0   
            end
        end
        
    else
        fill!(dc, zero(Float64))
        
    end
    
    @. concentrations = concentrations + dc;
    
end  