# 基于单一种群的特征选择方法
import warnings
warnings.filterwarnings('ignore')
import tensorflow as tf
import numpy as np
import pyes
import pandas as pd
import math
import matplotlib.pyplot as plt
from tensorflow.python.framework.errors_impl import InvalidArgumentError
from tensorflow.contrib.distributions import MultivariateNormalFullCovariance
from loguru import logger
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

tf.logging.set_verbosity(tf.logging.ERROR)

data = np.genfromtxt('../exampleData/glass.data', delimiter=',')

def encode(solution):
    """用于将个体编码成二进制形式，用于选择特征"""
    v_func = lambda x: 1.0 / (1.0 + np.exp(-x*2))
    return [True if v_func(x) > 0.5 else False for x in solution]

def DR(encoded_x):
    """计算特征维度缩减率"""
    return 1 - 1.0 * sum([1 for v in encoded_x if v]) / len(encoded_x)

def CA(encode_x, data_X, data_y):
    """计算分类准确率"""
    clf = make_pipeline(StandardScaler(), SVC(gamma='auto'))
    clf.fit(data_X[:,encode_x], data_y)
    return clf.score(data_X[:,encode_x], data_y)

# 目标函数包含 CA 和 DR 两部分内容，如果需要修改分类器，请修改 CA 函数
# 说明：文章中"图 4.2" 曲线与给定的公式 4.2 是对不上的，所以按最曲线重新定义了个 rho 函数
def rho(l):
    return 0.90 + 0.099 / (l/2000 + 1)

def objective(encode_x, data_X, data_y):
    l = len(encode_x)
    return rho(l) * CA(encode_x, data_X, data_y) + (1 - rho(l)) * DR(encode_x)

def real_objectives(lam, population, data_X, data_y):
    """
    适应值塑造,对应文章中的公式 4.3, 返回 lam 个塑造后的较优适应值，及各个体下标
    """
    encoded_population = [encode(x) for x in population]
    fitness_values = [objective(x, data_X, data_y) for x in encoded_population]
    # 取前 lam 个适应值较优的个体下标,再按升序排列
    selected = np.argsort(fitness_values)[::-1][:lam][::-1]

    temp = 1.0 * (fitness_values[selected[lam-1]] - fitness_values[selected[0]]) / (lam - 1) 
    
    real_fitness_values = [0] * lam

    for i, idx in enumerate(selected):
        if i == 0:
            real_fitness_values[i] = - (fitness_values[selected[lam-1]] - fitness_values[selected[i]]) / 2
        else:
            real_fitness_values[i] = real_fitness_values[0] + i * temp
    
    return real_fitness_values, selected

def event_on_generation(gen, population, current_best, history_best):
    best = history_best
    if gen == 0:
        best = current_best
    print('iter: {gen}, best solution(bool): {best}, best idxs: {idxs}'.format(gen=gen, best=encode(best), idxs='-'.join([str(i) for i, v in enumerate(encode(best)) if v])))

def NES(
    n_dimension: int, 
    objective, 
    n_iter, 
    learn_rate, 
    lam,
    mean,
    cov,
    random_state=None, 
    event_on_genration=event_on_generation):

    logger.debug({
        'dimension': n_dimension,
        'n_iter': n_iter,
        'learn_rate': learn_rate,
    })

    tf.set_random_seed(random_state)
    best, best_fitness = None, 1e+10
    
    def get_fitness(population): 
        return [-objective(solution) for solution in population]

    mean = tf.Variable(mean, dtype=tf.float32)
    cov = tf.Variable(cov, dtype=tf.float32)
    mvn = MultivariateNormalFullCovariance(loc=mean, covariance_matrix=cov)
    make_population = mvn.sample(lam)

    fitness_input = tf.placeholder(tf.float32, [lam, ])
    prob_output = tf.placeholder(tf.float32, [lam, n_dimension])
    loss = -tf.reduce_mean(mvn.log_prob(prob_output) * fitness_input)
    train_op = tf.train.GradientDescentOptimizer(learn_rate).minimize(loss)

    sess = tf.Session()
    sess.run(tf.global_variables_initializer())

    for g in range(n_iter):
        population = sess.run(make_population)

        scores = [objective(c) for c in population]
        
        event_on_genration(g, population, population[np.argsort(scores)[0]], best)

        fitness_values = get_fitness(population)
        sess.run(train_op, {fitness_input: fitness_values, prob_output: population}) 

        for i, score in enumerate([objective(x) for x in population]):
            if score < best_fitness:
                best, best_fitness = population[i], score

    return [best, best_fitness]

n_dimension = n_col - 1
n_iter = 10
learn_rate = 0.002
mu = 4
lam = 12
mean = data[:,:n_dimension].mean(axis=0)
cov = 10000 * tf.eye(n_dimension)
random_state = 1 # 固定输出

def es_objective(x):
    return -objective(encode(x), data[:,:n_dimension], data[:,n_dimension])

# 结论
# 原先 10 个特征，经过基于 NES 的特征选择方法，特征数缩减至 5 个，同时准确率保持在 95% 以上