<repo_name>GraceYYJ/cbx-k<file_sep>/ddpg_b.py
# encoding: utf-8
import numpy as np
import math
import my_common as mc

import torch
from quantization_env import QuantizationEnv
from ddpg import DDPG
from ounoise import OUNoise
from param_noise import AdaptiveParamNoiseSpec, ddpg_distance_metric
from replay_memory import ReplayMemory, Transition


class DDPGB(object):
    # x是数值向量
    # b是码值向量
    # c是标准codebook矩阵
    # action_output_num是码值输出的维度
    # replay_size是meomery队列的最大长度
    # new_b代表新计算得到的b
    # env代表action和obersevation的产生环境
    # agent代表实际的ddpg执行体
    # 保留这些噪声参数只是为了能够进入到需要随机探索的部分
    def __init__(self, C, b, x, action_output_num, actor_size, replay_size=1000000, ou_noise=True, param_noise=True,
                 noise_scale=0.3, final_noise_scale=0.3):
        self.C = C
        self.b = b
        self.x = x
        self.hd = action_output_num
        self.actor_size = actor_size
        self.memory = ReplayMemory(replay_size)
        self.new_b = None
        self.env = None
        self.agent = None
        self.ou_noise = ou_noise
        self.noise_scale = noise_scale
        self.final_noise_scale = final_noise_scale
        self.ounoise = OUNoise(action_output_num) if ou_noise else None
        self.param_noise = AdaptiveParamNoiseSpec(initial_stddev=0.05, desired_action_stddev=noise_scale,
                                                  adaptation_coefficient=1.05) if param_noise else None

    def update_B(self, c, b, x):
        self.C = c
        self.b = b
        self.x = x

    # 备选coff代表reward中的权重比例[0.2, 0.8]
    def generate_B(self, coff, gamma, tau, hidden_size, num_inputs, actor_size, num_episodes=60000,
                   exploration_end=150, batch_size=512, updates_per_step=5000):

        self.env = QuantizationEnv(self.C, self.b, self.x, self.hd, coff)
        self.agent = DDPG(gamma, tau, hidden_size, self.env.action_bin, num_inputs, actor_size)
        rewards = []
        total_numsteps = 0
        updates = 0
        max_trail = 10000
        best_bb = 10000

        # 开启num_episodes次最佳方案寻找
        for i_episode in range(num_episodes):
            state = torch.Tensor([self.env.reset()])
            if self.ou_noise:
                self.ounoise.scale = (self.noise_scale - self.final_noise_scale) * max(0, exploration_end - i_episode) \
                                     / exploration_end + self.final_noise_scale
                self.ounoise.reset()
            if self.param_noise:
                self.agent.perturb_actor_parameters(self.param_noise)

            episode_reward = 0
            continuous_neg = 0
            continuous_pos = 0
            temp_trail = 0

            control_bit = 0
            next_state = self.env.compute_Cbx(self.b)
            next_state = torch.Tensor([next_state])

            while True:
                # yyj
                if control_bit > 15:
                    control_bit = control_bit % 16
                state = next_state
                action = self.agent.select_action(state, self.ounoise, self.param_noise)
                next_state, reward, done, bb = self.env.step(action, control_bit,self.actor_size)
                # print(control_bit, next_state[0], reward, done, bb)
                control_bit = control_bit + 1
                total_numsteps += 1
                episode_reward += reward
                # bb是c_v值
                if best_bb > bb:
                    best_bb = bb
                    self.new_b = action

                if reward > 0:
                    continuous_pos += 1
                    continuous_neg = 0
                    if continuous_pos > 10:
                        done = True
                if reward < 0:
                    continuous_neg += 1
                    continuous_pos = 0
                    if continuous_neg > 10:
                        done = True
                if temp_trail > max_trail:
                    done = True

                action = torch.Tensor(action)
                mask = torch.Tensor([not done])
                next_state = torch.Tensor([next_state])
                reward = torch.Tensor([reward])

                self.memory.push(state, action, mask, next_state, reward)

                # state = next_state
                temp_trail += 1

                # memorysize还不够，不会进入
                if len(self.memory) > batch_size:
                    for _ in range(updates_per_step):
                        transitions = self.memory.sample(1)
                        batch = Transition(*zip(*transitions))

                        # value_loss属于右边的网络，policy_loss属于左边的网络
                        value_loss, policy_loss = self.agent.update_parameters(batch)
                        print("epoch:", i_episode, "updates", updates, "value_loss:",
                              value_loss, " policy_loss:", policy_loss)

                        updates += 1
                if done:
                    break

            if self.param_noise:
                episode_transitions = self.memory.memory[self.memory.position - batch_size:self.memory.position]
                states = torch.cat([transition[0] for transition in episode_transitions], 0)
                unperturbed_actions = self.agent.select_action(states, None, None)
                perturbed_actions = torch.cat([transition[1] for transition in episode_transitions], 0)

                ddpg_dist = ddpg_distance_metric(perturbed_actions.numpy(), unperturbed_actions.numpy())
                self.param_noise.adapt(ddpg_dist)

            rewards.append(episode_reward)
            continuous_neg = 0
            continuous_pos = 0
            temp_trail = 0
            if i_episode % 10 == 0 and i_episode != 0:
                state = torch.Tensor([self.env.reset()])
                episode_reward = 0
                control_bit = 0
                while True:
                    action = self.agent.select_action(state)
                    next_state, reward, done, bb = self.env.step(action.numpy()[0], control_bit)
                    episode_reward += reward
                    if best_bb > bb:
                        best_bb = bb
                        self.new_b = action

                    if reward > 0:
                        continuous_pos += 1
                        continuous_neg = 0
                        if continuous_pos > 10:
                            done = True
                    if reward < 0:
                        continuous_neg += 1
                        continuous_pos = 0
                        if continuous_neg > 10:
                            done = True
                    if temp_trail > max_trail:
                        done = True

                    next_state = torch.Tensor([next_state])

                    state = next_state
                    temp_trail += 1

                    if done:
                        break

                rewards.append(episode_reward)
                print(
                    "Episode: {}, total numsteps: {}, reward: {}, average reward: {}".format(i_episode, total_numsteps,
                                                                                             rewards[-1],
                                                                                             np.mean(rewards[-10:])))

        return self.new_b
<file_sep>/Functions.py
#!/usr/bin/env python3
# coding=utf-8

from __future__ import absolute_import
from __future__ import print_function, division
import os
from os import path
import numpy as np
import sys
import random
from datetime import datetime

def PrintWithTime(string=None):
    if string is None:
        print("{0}".format(datetime.now().strftime("%m-%d %H:%M:%S")))
    else:
        print("{0}  {1}".format(datetime.now().strftime("%m-%d %H:%M:%S"), string))


def ProgressBar(progress):
    progress *= 100
    progress = round(progress)
    a = '=' * int(progress) + '-' * (100 - int(progress))
    a = a[:47] + (' 100 ' if progress >= 100 else " %2d%% " %
                                                  progress) + a[52:] + '\n'
    print(a)


def BarFormat(string):
    l = len(string)
    a = (100 - 2 - l) // 2
    b = 2 * a + 2 + l == 99
    return '=' * a + ' {0} '.format(string) + '=' * (a + 1) if b else '=' * a + ' {0} '.format(string) + '=' * a

def ObjectFunction(C, b, x):
    # || Cb - x ||
    # [K, ] * [K, D] - [D, ]
    return np.linalg.norm(b.dot(C) - x)


def InitializeBinary(k):
    """Init a random b (0, 1) ^ [K, 1] for first iteration

    Args:
        k (int): the num of b
    """
    b = np.zeros([k, ], dtype=np.int8)

    """ Primal test, half of b set to 1, and another half is 0 """
    choice = np.random.choice(k, k // 2)
    b[choice] = 1
    return b


def OptimizeC(C, X, B, learning_rate=1e-2):
    """[summary]
    
    Args:
        C (ndarray): Codebook, [K, D]
        X (ndarray): Batch of X, [N, D]
        B (ndarray): Code of this batch, [N, K]
    """
    delta = 2 * B.T.dot(B.dot(C) - X)
    return C - learning_rate * delta



def MiniBatchOptimizeC(C, X, B, learning_rate=1e-2):
    """[summary]
    
    Args:
        C (ndarray): Codebook, [K, D]
        X (ndarray): Batch of X, [N, D]
        B (ndarray): Code of this batch, [N, K]
    """
    N = X.shape[0]
    delta = 2 * B.T.dot(B.dot(C) - X)
    return C - learning_rate * delta / N<file_sep>/test.py
#encoding: utf-8
import numpy as np
import random
'''
for Age_laoshi in range(5,100):
    for Age_gap in range(5,100):
        if Age_laoshi+Age_gap+Age_laoshi == 51 and Age_laoshi- Age_gap == 4 * (Age_laoshi- Age_gap - Age_gap):
            print Age_laoshi
            print Age_gap
            break


for a in range(1,10):
    for b in range(0,10):
        for c in range(0,10):
            x = 100*a + 10*b + c
            if x//11 - x % 11 == 1:
                print x

for black in range(1,1000):
    for white in range(1,1000):
        if 5*black+white == 2*(black+white):
            print (5*white+black)/(black+white)
'''


# to_tensor = ToTensor()
# to_pil = ToPILImage()
# #
# lena = Image.open("imgs/lena.png")
# input = to_tensor(lena).unsqueeze(0)
# # print(input.size())
# # kernel = t.ones(3, 3, 3) / -9
# # kernel[1][1] = 1
# # conv = nn.Conv2d(1, 1, (3, 3), 1, bias=False)
# # conv.weight.data = kernel.view(1, 3, 3, 3)
# # out = conv(V(input))
# # print(to_pil(out.data.squeeze(0)).show())
#
# pool = nn.AvgPool2d(2,2)
# out = pool(V(input))
# print(out.size())
# print(to_pil(out.data.squeeze(0)).show())


#input = V(t.randn(2, 3))
#linear = nn.Linear(3, 4)
#h = linear(input)
# print(h.size())
# print(h)
# print(sum(h))

# bn = nn.BatchNorm1d(4) #批标准化，4通道，(mear,std):(0,1)->(0,4)
# bn.weight.data = t.ones(4) * 4 #设置参数
# bn.bias.data = t.zeros(4) #设置参数
# bn_out = bn(h)
# # print(bn_out.mean(0))
# # print(bn_out.var(0, unbiased=True))
# dropout = nn.Dropout(0.8) #对输入的每个元素以0.8的概率置为0
# o = dropout(bn_out)
# print(o)

#a = np.array([i % 2 for i in np.random.randint(1, 100, size=16)])
#print a

list = random.sample(range(0, 15), 8)
vec = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
for i in list:
    vec[i] = -vec[i]
print vec

a = False
b = 5
print a*b<file_sep>/tabusearch.py
#!/usr/bin/env python3
# coding=utf-8

from __future__ import absolute_import
from __future__ import print_function, division
import os
from os import path
import numpy as np
# from sift_read import sift10k_read
from ITS import ITS_Instance
from Functions import *
from multiprocessing import Pool
from time import perf_counter

# 在这里用的是BC-X，我们用的是CB-X
# K = 256
K = 16
# D = 64
D = 128
batch_size = 1
Total = batch_size * 1
jump_count = 4

Codes = np.zeros([Total, K], dtype=np.int)

# [k, D]
C = np.random.rand(K, D)

b = np.random.rand(Total, K)

bits = np.zeros([Total, K], dtype=np.int)

bits[b > 0.5] = 1

X = bits.dot(C)

X = X + np.random.randn(Total, D)


def partition(arg):
    index, Query = arg
    distortion = 0
    N = Query.shape[0]
    codes = np.zeros([N, K], dtype=np.int)
    i = 0
    for q in Query:
        tabu = ITS_Instance(5 * D, jump_count, C, q)
        tabu.run()
        codes[i] = tabu.BestSolution
        d = tabu._fit(tabu.BestSolution)
        distortion += d
        i += 1
        # print((index, i))
    print(index, 'Completed')
    return distortion, codes


learning_rate = 1e-2

curve = []

def tabuSearch(tabu_t):
    global C,learning_rate,curve
    for _ in range(tabu_t):
        distortion = 0

        np.random.shuffle(X) #将X按行重新排列

        splitted = np.split(X, Total // batch_size) ##将X分成Total//batch_size个组

        l = list()
        i = 0
        for s in splitted:
            l.append((i, s))
            i += 1
        with Pool(4) as p:
            results = p.map(partition, l) #l为list，对l中每个元素进行partition运算，一次运算四个元素，因为只有4个进程
            k = 0
            i = 0
            for res in results:
                dis, code = res
                n = code.shape[0]
                distortion += dis
                Codes[k:k+n] = code #按组更新
                C = MiniBatchOptimizeC(C, splitted[i], code, learning_rate)
                i += 1
                k += n
        distortion = distortion / X.shape[0]
        PrintWithTime(distortion)
        curve += distortion
        learning_rate *= 0.99
    # print("C:", C)
    # print("C.shape:", C.shape)
    # print("Codes:", Codes)
    # print("Codes.shape:", Codes.shape)
    # print("curve:", curve)
    # print("X.shape:",X.shape)
    # curve = np.array(curve)
    # print("np.array(curve):", np.array(curve))
    # exit()
    # np.save('./distortion_curve', curve)
    return C.T,Codes.T
<file_sep>/ddpg.py
# encoding: utf-8
import sys

import torch
import math
import torch.nn as nn
from torch.optim import Adam
from torch.autograd import Variable
import torch.nn.functional as F
import my_common as mc
from spp_layer import spatial_pyramid_pool


def soft_update(target, source, tau):
    for target_param, param in zip(target.parameters(), source.parameters()):
        target_param.data.copy_(target_param.data * (1.0 - tau) + param.data * tau)


def hard_update(target, source):
    for target_param, param in zip(target.parameters(), source.parameters()):
        target_param.data.copy_(param.data)


spp_num_outputs = [1, 2, 4]


class ActorRnn(nn.Module):
    def __init__(self, in_dim, hidden_dim, n_layer, n_class):
        super(ActorRnn, self).__init__()
        self.n_layer = n_layer
        self.hidden_dim = hidden_dim
        self.lstm = nn.LSTM(in_dim=in_dim, hidden_dim=hidden_dim, n_layer=n_layer, batch_first=True)

    def forward(self, C):
        out, _ = self.lstm(C)
        # res_Cbx = states
        # res_Cbx = res_Cbx.squeeze()

        # Decode hidden state of last time step
        out = self.fc(out[:, -1, :])

        return out


class Actor(nn.Module):
    # hidden_size代表全连接中的隐层基本数量
    # action代表b
    # num_outputs代表Actor网络的输出维度
    # num_intputs代表Actor网络的输入维度

    # -------------------- 如果使用resnet 则使用下列参数 -------------------------
    # spp_num_outputs代表数据使用spp时的分割字典数组
    # spp_data_width = data_width代表数据从一位向量被转化为二维向量时的第一个维度数
    # fc_num_inputs代表在spp下计算得到的金字塔输出维度，也是fc的输入维度，取代num_intputs
    def __init__(self, hidden_size, action, num_inputs, num_output, spp_num_outputs=[1, 2, 4], data_width=8):
        super(Actor, self).__init__()
        self.action = action
        self.num_inputs = num_inputs
        self.num_outputs = num_output

        self.spp_data_width = data_width
        self.spp_num_outputs = spp_num_outputs
        self.fc_num_inputs = sum([i * i for i in self.spp_num_outputs])

        '''
        #设计残差模块用于Cb-x的向量输入
        self.Cbx_res_block1 = nn.Sequential(
            nn.Conv2d(1, 128, kernel_size=1, padding=0),
            nn.BatchNorm2d(128),
            nn.RReLU(0.66, 0.99),
            nn.Conv2d(128, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.RReLU(0.33, 0.66),
            nn.Conv2d(128, 1, kernel_size=3, padding=1),
            nn.BatchNorm2d(1)
        )

        self.Cbx_res_block2 = nn.Sequential(
            nn.Conv2d(1, 128, kernel_size=1, padding=0),
            nn.BatchNorm2d(128),
            nn.RReLU(0.66, 0.99),
            nn.Conv2d(128, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.RReLU(0.33, 0.66),
            nn.Conv2d(128, 1, kernel_size=3, padding=1),
            nn.BatchNorm2d(1)
        )
        '''

        # 第一个全连接网络，层上归一化，随机relu
        # self.afc1 = nn.Linear(self.fc_num_inputs, hidden_size)
        self.afc1 = nn.Linear(self.num_inputs, hidden_size)
        self.ln1 = nn.LayerNorm(hidden_size)
        self.relu1 = nn.RReLU(0.01, 0.33)

        # 第二个全连接网络，层上归一化，随机relu
        self.afc2 = nn.Linear(hidden_size, hidden_size)
        self.ln2 = nn.LayerNorm(hidden_size)
        self.relu2 = nn.RReLU(0.33, 0.66)

        # 第三个全连接网络，层上归一化，随机relu
        self.afc3 = nn.Linear(hidden_size, hidden_size)
        self.ln3 = nn.LayerNorm(hidden_size)
        self.relu3 = nn.RReLU(0.66, 0.99)

        # 第四个全连接网络，层上归一化，随机relu
        self.afc4 = nn.Linear(hidden_size, hidden_size // 2)
        self.ln4 = nn.LayerNorm(hidden_size // 2)
        self.drop = nn.Dropout()

        # 第五层全连接，准备进行激活
        self.mu = nn.Linear(hidden_size // 2, self.num_outputs)
        self.mu.weight.data.mul_(0.1)
        self.mu.bias.data.mul_(0.1)

    def forward(self, states):
        res_Cbx = states
        '''
        res_Cbx = res_Cbx.unsqueeze(2)
        res_Cbx = res_Cbx.unsqueeze(3)
        res_Cbx = res_Cbx.view(res_Cbx.size(0), 1, self.spp_data_width, res_Cbx.size(1) // self.spp_data_width)

        x_Cbx = self.Cbx_res_block1(res_Cbx)
        x_Cbx += res_Cbx
        res_Cbx = x_Cbx
        x_Cbx = self.Cbx_res_block2(res_Cbx)
        x_Cbx += res_Cbx

        x_Cbx = spatial_pyramid_pool(x_Cbx, x_Cbx.size(0), [x_Cbx.size(2), x_Cbx.size(3)], self.spp_num_outputs)
        x = x_Cbx.squeeze(1)
        '''
        res_Cbx = res_Cbx.squeeze()

        x = self.afc1(res_Cbx)
        x = self.ln1(x)
        x = self.relu1(x)

        x = self.afc2(x)
        x = self.ln2(x)
        x = self.relu2(x)

        x = self.afc3(x)
        x = self.ln3(x)
        x = self.relu3(x)

        x = self.afc4(x)
        x = self.ln4(x)
        x = self.drop(x)

        mu = torch.sigmoid(self.mu(x))
        return mu


class Critic(nn.Module):
    # 与Actor网络几乎一致，只是在前两层的网络中
    def __init__(self, hidden_size, action, num_inputs, spp_num_outputs, data_width=8):
        super(Critic, self).__init__()
        self.action = action
        self.num_outputs = self.action.shape[0]
        # self.num_outputs=1
        self.num_inputs = num_inputs

        self.spp_data_width = data_width
        self.spp_num_outputs = spp_num_outputs
        self.fc_num_inputs = sum([i * i for i in self.spp_num_outputs])

        '''
        self.Cbx_res_block1 = nn.Sequential(
            nn.Conv2d(1, 128, kernel_size=1, padding=0),
            nn.BatchNorm2d(128),
            nn.RReLU(0.66, 0.99),
            nn.Conv2d(128, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.RReLU(0.33, 0.66),
            nn.Conv2d(128, 1, kernel_size=3, padding=1),
            nn.BatchNorm2d(1)
        )

        self.Cbx_res_block2 = nn.Sequential(
            nn.Conv2d(1, 128, kernel_size=1, padding=0),
            nn.BatchNorm2d(128),
            nn.RReLU(0.66, 0.99),
            nn.Conv2d(128, 128, kernel_size=3, padding=1),
            nn.BatchNorm2d(128),
            nn.RReLU(0.33, 0.66),
            nn.Conv2d(128, 1, kernel_size=3, padding=1),
            nn.BatchNorm2d(1)
        )
        '''

        self.cfc11 = nn.Linear(self.num_outputs, hidden_size)
        self.ln11 = nn.LayerNorm(hidden_size)
        self.relu11 = nn.RReLU(0.01, 0.33)

        self.cfc12 = nn.Linear(self.num_inputs, hidden_size)
        self.ln12 = nn.LayerNorm(hidden_size)
        self.relu12 = nn.RReLU(0.01, 0.33)

        self.cfc2 = nn.Linear(hidden_size * 2, hidden_size * 2)
        self.ln2 = nn.LayerNorm(hidden_size * 2)
        self.relu2 = nn.RReLU(0.33, 0.66)

        self.cfc3 = nn.Linear(hidden_size * 2, hidden_size)
        self.ln3 = nn.LayerNorm(hidden_size)
        self.relu3 = nn.RReLU(0.66, 0.99)

        self.cfc4 = nn.Linear(hidden_size, hidden_size // 2)
        self.ln4 = nn.LayerNorm(hidden_size // 2)
        self.drop = nn.Dropout()

        self.V = nn.Linear(hidden_size // 2, 1)
        self.V.weight.data.mul_(0.1)
        self.V.bias.data.mul_(0.1)

    def forward(self, states, actions):
        states = states.squeeze()
        actions = actions.squeeze()

        x_b = actions
        x_b = x_b.reshape(-1, self.num_outputs)

        x_b = self.cfc11(x_b)
        x_b = self.ln11(x_b)
        x_b = self.relu11(x_b)

        x_Cbx = states.squeeze()
        x_Cbx = self.cfc12(x_Cbx)
        x_Cbx = self.ln12(x_Cbx)
        x_Cbx = self.relu12(x_Cbx)

        x = torch.cat([x_b, x_Cbx], -1)
        # x = x.squeeze(1)

        x = self.cfc2(x)
        x = self.ln2(x)
        x = self.relu2(x)

        x = self.cfc3(x)
        x = self.ln3(x)
        x = self.relu3(x)

        x = self.cfc4(x)
        x = self.ln4(x)
        x = self.drop(x)

        V = self.V(x)
        return V


class DDPG(object):
    def __init__(self, gamma, tau, hidden_size, action, num_inputs, actor_size):

        self.action = action
        self.hd = action.shape[0]

        if torch.cuda.is_available():
            self.actor = Actor(hidden_size, self.action, num_inputs, actor_size).cuda()
            self.actor_target = Actor(hidden_size, self.action, num_inputs, actor_size).cuda()
            self.actor_perturbed = Actor(hidden_size, self.action, num_inputs, actor_size).cuda()  # 探索阶段
        else:
            self.actor = Actor(hidden_size, self.action, num_inputs, actor_size)
            self.actor_target = Actor(hidden_size, self.action, num_inputs, actor_size)
            self.actor_perturbed = Actor(hidden_size, self.action, num_inputs, actor_size)  # 探索阶段

        self.actor_optim = Adam(self.actor.parameters(), lr=0.05)  # actor的学习阶段，随机梯度下降，注意学习率

        if torch.cuda.is_available():
            self.critic = Critic(hidden_size, self.action, num_inputs, spp_num_outputs).cuda()
            self.critic_target = Critic(hidden_size, self.action, num_inputs, spp_num_outputs).cuda()
        else:
            self.critic = Critic(hidden_size, self.action, num_inputs, spp_num_outputs)
            self.critic_target = Critic(hidden_size, self.action, num_inputs, spp_num_outputs)

        self.critic_optim = Adam(self.critic.parameters(), lr=0.1)  # critic的学习阶段，随机梯度下降，注意学习率

        self.gamma = gamma
        self.tau = tau

        hard_update(self.actor_target, self.actor)  # Make sure target is with the same weight
        hard_update(self.critic_target, self.critic)

    # 探索阶段的行为确定，对mu的扰动
    def select_action(self, state, action_noise=None, param_noise=None):
        self.actor.eval()
        if param_noise is not None:
            if (torch.cuda.is_available()):
                mu = self.actor_perturbed((Variable(state).cuda()))
            else:
                mu = self.actor_perturbed((Variable(state)))
        else:
            if (torch.cuda.is_available()):
                mu = self.actor((Variable(state).cuda()))
            else:
                mu = self.actor((Variable(state)))

        self.actor.train()
        mu = mu.data

        # 扰动的出发条件
        if action_noise is not None:
            '''
            if torch.cuda.is_available():
                mu += torch.Tensor(action_noise.noise()).cuda()
            else:
                mu += torch.Tensor(action_noise.noise())
            '''
            # 根据随机的效果取一半的数据更改符号
            mu = mc.random_change_sgn(0, self.hd - 1, self.hd // 22, mu)

        return mu

    def update_parameters(self, batch):
        if torch.cuda.is_available():
            state_batch = Variable(torch.cat(batch.state)).cuda()
            action_batch = Variable(torch.cat(batch.action)).cuda()
            reward_batch = Variable(torch.cat(batch.reward)).cuda()
            mask_batch = Variable(torch.cat(batch.mask)).cuda()
            next_state_batch = Variable(torch.cat(batch.next_state)).cuda()
        else:
            state_batch = Variable(torch.cat(batch.state))
            action_batch = Variable(torch.cat(batch.action))
            reward_batch = Variable(torch.cat(batch.reward))
            mask_batch = Variable(torch.cat(batch.mask))
            next_state_batch = Variable(torch.cat(batch.next_state))

        next_action_batch = self.actor_target(next_state_batch)
        next_state_action_values = self.critic_target(next_state_batch, next_action_batch)

        reward_batch = reward_batch.unsqueeze(1)  # 纵向排列
        mask_batch = mask_batch.unsqueeze(1)
        expected_state_action_batch = reward_batch + (
                self.gamma * mask_batch * next_state_action_values)  # mask_batch中是true和false，python会自动将其转换为1和0

        # critic网络重置梯度为0
        self.critic_optim.zero_grad()

        state_action_batch = self.critic((state_batch), (action_batch))

        # 以均方差loss计算value值并反馈
        value_loss = F.mse_loss(state_action_batch, expected_state_action_batch)
        value_loss.backward()
        self.critic_optim.step()  # 完成参数更新

        # critic网络重置梯度为0
        self.actor_optim.zero_grad()

        # 对于Actor网络，使用critic的计算值作为反馈的loss
        policy_loss = -self.critic((state_batch), self.actor((state_batch)))
        policy_loss = policy_loss.mean()
        policy_loss.backward()
        self.actor_optim.step()

        soft_update(self.actor_target, self.actor, self.tau)
        soft_update(self.critic_target, self.critic, self.tau)

        return value_loss.item(), policy_loss.item()

    def perturb_actor_parameters(self, param_noise):
        """Apply parameter noise to actor model, for exploration"""
        hard_update(self.actor_perturbed, self.actor)
        params = self.actor_perturbed.state_dict()  # 取出所有的params，并且把给的noise加入到param中
        for name in params:
            if 'ln' in name:
                pass
            param = params[name]
            param += torch.randn(param.shape) * param_noise.current_stddev

    def save_model(self, env_name, suffix="", actor_path=None, critic_path=None):
        if not os.path.exists('models/'):
            os.makedirs('models/')

        if actor_path is None:
            actor_path = "models/ddpg_actor_{}_{}".format(env_name, suffix)
        if critic_path is None:
            critic_path = "models/ddpg_critic_{}_{}".format(env_name, suffix)
        print('Saving models to {} and {}'.format(actor_path, critic_path))
        torch.save(self.actor.state_dict(), actor_path)
        torch.save(self.critic.state_dict(), critic_path)

    def load_model(self, actor_path, critic_path):
        print('Loading models from {} and {}'.format(actor_path, critic_path))
        if actor_path is not None:
            self.actor.load_state_dict(torch.load(actor_path))
        if critic_path is not None:
            self.critic.load_state_dict(torch.load(critic_path))
<file_sep>/quantization_env.py
# encoding: utf-8
import numpy as np
import my_common as mc
from torch.autograd import Variable
import math
import torch

class QuantizationEnv(object):
    #codebook是np格式的矩阵m*n
    #b以01向量进行计算的结果为hard
    #b以浮点数向量进行计算的结果为soft
    #self.action一定是01向量
    def __init__(self, codebook, b, x, code_len, coff):
        self.reward_coff = coff
        self.b = b
        self.x = x
        self.C = codebook
        self.cl = code_len
        #yyj
        self.action_bin=self.b

        self.observation = self.generate_state(self.action_bin)
        self.mask = False
        self.hard_pre_v = 0
        self.hard_c_v = 0
        self.soft_pre_v = 0
        self.soft_c_v = 0
        self.is_first_v = True
        self.reward = 0
        self.history_list = []
        self.info = []

    #还原状态为初始状态
    def reset(self):
        #yyj
        self.action_bin = self.b
        self.observation = self.generate_state(self.action_bin)
        self.mask = False
        self.hard_pre_v = 0
        self.hard_c_v = 0
        self.soft_pre_v = 0
        self.soft_c_v = 0
        self.is_first_v = True
        self.reward = 0

        return self.observation

    #step by yyj
    def step(self, action, control_bit, actor_size):
        #如果得到的b是全0的，则需要重新随机一个b作为action

        self.action_bin[control_bit:control_bit+actor_size] = mc.binazation(action)
        self.action_b[control_bit:control_bit+actor_size] = action.numpy()
        #action_b=self.b
        #action_b[control_bit]=action
        if self.action_bin.sum() == 0:
            #self.action = self.rand_action(self.cl)
            #全0可以看做是出现了一种错误的情况
            self.mask = True

        self.observation = self.compute_Cbx(self.action_bin)
        self.soft_observation=self.compute_Cbx(self.action_b)

        if self.is_first_v:
            #yyj 第一次的时候 pre是用随机b所计算出的cb-x的模
            self.hard_pre_v=self.compute_v(self.compute_Cbx(self.b))
            self.hard_c_v = self.compute_v(self.observation)
            self.soft_pre_v = self.compute_v(self.compute_Cbx(self.b))
            self.soft_c_v = self.compute_v(self.soft_observation)
            self.is_first_v = False
        else:
            self.hard_pre_v = self.hard_c_v
            self.hard_c_v = self.compute_v(self.observation)
            self.soft_pre_v = self.soft_c_v
            self.soft_c_v = self.compute_v(self.soft_observation)
            self.reward = self.compute_reward()

        if self.hard_c_v < 10 or self.hard_c_v > 10000:
            self.mask = True
        #如果两次计算得到的浮点数向量很靠近，说明已经收敛可以跳出了
        if (self.soft_pre_v - self.soft_c_v).sum() < 0.0001:
             self.mask = True
        return self.observation, self.reward, self.mask, self.hard_c_v

    #计算reward
    # 如果当前得到的模长大于上一次得到的模长，说明这个b的走势不好，reward应当给负值
    # 反之，应当是正值
    def compute_reward(self):
        #np.linalg.norm(a,2)所有元素的平方和开更号
        hard = (self.hard_pre_v - self.hard_c_v)/self.hard_pre_v
        soft = (self.soft_pre_v - self.soft_c_v)/self.soft_pre_v
        return self.reward_coff[0] * hard + self.reward_coff[1] * soft

    #计算矩阵C与b的乘积
    def compute_Cb(self, b):
        return np.dot(self.C, b)

    #计算Cb-x
    def compute_Cbx(self, b):
        return self.compute_Cb(b) - self.x

    #统计得到当前cb-x的模长（这个值采用绝对值加和，如果越小，说明b的效果越好）
    def compute_v(self, observation):
        Cbx = observation
        sum_Cbx = np.linalg.norm(Cbx, 1)
        return sum_Cbx

    #随机一个cl长度的01向量
    def rand_action(self, cl):
        return mc.random_vecb(cl)

    # 将Cb-x作为状态空间
    def generate_state(self, b):
        return self.compute_Cbx(b)
<file_sep>/main.py
# encoding: utf-8
import numpy as np
import random_c as rc
import my_common as mc
import gradient_c as gc
import ddpg_b as db
import math
import torch
import fvecs_read
from tabusearch import *
#问题综述：找到X的最佳量化模式，数学表达为求解CB-X的最小值，C和B迭代进行，C为根据维度的随机初始矩阵
#X是数据组成的矩阵，维度为d*n
#C是codebook，是一个d*k的矩阵
#B是编码矩阵,并且元素全为0和1，维度k*n

#解决思路：
# 1、在使用禁忌搜索的方法迭代z(z=5)轮之后，得到计算过程中相对靠谱的C
# 2、使用强化学习方案，随机选取X中的任意t(t=5)个向量作为学习样本,使用较大的epoach次数与step次数，完成t个样本的最佳b路径选择模型
# 3、对剩余的n-t个样本，依次对比t个样本，选取一范式差值最小的t所在的模型，进行finetune运算，得到各自的b值
# 4、整理n个b值形成B后，使用梯度下降（矩阵求导）的方式得到新的C
# 5、返回2进行一下轮迭代，直到迭代次数终止
# 6、收集最后一次t个样本及模型，用于可能到来的查询

#参数设定
static_d = 0
static_n = 0
static_k = 16
org_X = []
static_X = []
X_loss = 0
_C = []
_B = []
seed_c = 4
X_loss_thr = 88
#禁忌搜索的最大迭代次数
tabu_t = 1 #5
#_C与_B的循环的最大迭代次数
Max_iter_num = 18
#输出的b的维度
action_output_num = 16
#actor网络输出维度
actor_size=4
coff = [0.8, 0.2]
gamma = 0.99
tau = 0.001
hidden_size = 128

#读取所有的X
# org_X = np.loadtxt(open("mock_data", "rb"), delimiter=" ", skiprows=0)
org_X = fvecs_read.fvecs_read("./datasets/siftsmall/siftsmall_base.fvecs")
# print(org_X.shape)
static_X = org_X.T
static_n = org_X.shape[0]
static_d = static_X.shape[0]

#在禁忌搜索z轮迭代下进行强化学习的_C和_B得到
#进入强化学习的迭代循环
temp_iter_num = tabu_t
#从X中随机得到t个样本(以多进程的方式对他们进行并行的模型训练)

if __name__ == "__main__":
    # 使用禁忌搜索得到_C和_B
    _C, _B = tabuSearch(tabu_t)
    i = 0
    while True:
        #对抽取的t个样本中的第一个样本进行学习
        # _C是当前codebook，b是这个样本对应的码值(来源于上一阶段的禁忌搜索)，x是抽取到的样本
        # model_ddpg1 = db.DDPGB(_C, b, x, action_output_num)
        model_ddpg1 = db.DDPGB(_C, _B, org_X[i:i+1].T, action_output_num,actor_size)
        float_b1 = model_ddpg1.generate_B(coff, gamma, tau, hidden_size, static_d,actor_size)
        b1 = mc.binazation(float_b1)
        # 对于剩余的x样本，比对和这t个样本的欧式距离，选取距离最小的进行模型读取和finetune，得到各自的b
        # 组合这些b形成_B
        # 根据得到的B和X通过矩阵求导和梯度下降法得到新的C
        obj_gc = gc.GrandientC(_C, _B, static_X)
        _C, _C_loss = obj_gc.degrade(0.008, math.pow(10, -9))

        #检验结果
        X_loss = mc.matrix_distance(np.dot(_C, _B) - static_X)
        print("the {} epoach loss is {}".format(temp_iter_num, X_loss))

        temp_iter_num += 1
        if X_loss < X_loss_thr:
            break
        if temp_iter_num > Max_iter_num:
            break
        i+=1

    #给出最后的C和B
    print("_C:",_C)

    print("_B:",_B)
<file_sep>/gradient_c.py
# encoding: utf-8
import numpy as np
import my_common as mc

class GrandientC(object):
    def __init__(self, c, b, x):
        self.C = c
        self.B = b
        self.X = x
        self.loss = []
        print("degrade for C is ready")

    def degrade(self, lamda, epsion, maxiter=10000):
        _maxiter = maxiter
        print("degrade for C is start")
        while True:
            _maxiter = _maxiter - 1
            delta_C = self.deltaC()
            loss = mc.matrix_distance(lamda * delta_C)
            # print(loss)
            self.loss.append(loss)
            self.C = self.C - lamda * delta_C
            if loss < epsion or _maxiter == 0:
                return self.C, self.loss

    def deltaC(self):
        return 2 * np.dot((np.dot(self.C, self.B) - self.X), self.B.T)

<file_sep>/ITS.py
﻿#!/usr/bin/env python3
# coding=utf-8

from __future__ import absolute_import
from __future__ import print_function, division
import os
from os import path
import numpy as np
import sys
import random
from time import perf_counter


class ITS(object):
    def __init__(self, maxIteration: int, maxRunningTicks: int, initFunc, fitFunc):
        """禁忌搜索初始化

        Args:
            maxIteration (int): tabu search 迭代次数
            maxRunningTicks (int): 停机时间, 秒
            initFunc ([type]): 解的初始化函数
            fitFunc ([type]): 目标评价函数
        """

        self.BestSolution = None
        self.BestFitness = 0.0

        self._initSolution = initFunc
        self._fit = fitFunc

        self._maxIteration = maxIteration
        self._maxRunningTicks = maxRunningTicks

    def run(self):

        startTime = perf_counter()

        currentSolution = self._initSolution()

        D = currentSolution.shape[0]
        oneThird = D // 3

        currentFitness = self._fit(currentSolution)

        tabuList = np.zeros([D], dtype=np.float)

        self.BestFitness = currentFitness
        self.BestSolution = currentSolution

        # 非禁忌候选集中的解
        nextSolution = None
        nextFitness = int(sys.maxsize)
        nextTabu0 = None
        fitList = np.zeros([D], dtype=np.float)

        # 禁忌候选集中的解
        tabuSolution = None
        netxtTabuFitness = int(sys.maxsize)
        nextTabu1 = None
        tabufitList = np.zeros([D], dtype=np.float)

        total_iter = 0

        while True:
            bestList = None
            tabuList = np.zeros([D], dtype=np.float)
            # 迭代限制
            for i in range(self._maxIteration):
                nextFitness = int(sys.maxsize)
                netxtTabuFitness = int(sys.maxsize)

                success = False

                # 邻域候选集中查找
                for j in range(D):
                    neighbor = np.copy(currentSolution)
                    neighbor[j] = 1 - neighbor[j]
                    val = self._fit(neighbor)
                    # 是否满足禁忌要求
                    if i > tabuList[j]:
                        fitList[j] = val
                        # 选择邻域中评价最好的解作为候选
                        if nextFitness > val:
                            nextSolution = neighbor
                            nextFitness = val
                            # 记录位置
                            nextTabu0 = j
                            success = True
                    else:
                        tabufitList[j] = val
                        # 选择邻域中评价最好的解作为候选
                        if netxtTabuFitness > val:
                            tabuSolution = neighbor
                            netxtTabuFitness = val
                            # 记录位置
                            nextTabu1 = j

                # 赦免：禁忌表已满，或禁忌中的值比最好值还好
                if not success or (nextFitness > netxtTabuFitness and netxtTabuFitness < self.BestFitness):
                    # 更新禁忌表，当前值
                    # tabuList[nextTabu1] = i + (D / 150) + random.randint(1, 10)
                    tabuList = np.zeros([D], dtype=np.float)
                    currentSolution = tabuSolution
                    currentFitness = netxtTabuFitness
                    fitnessList = tabufitList
                else:
                    # 更新禁忌表，当前值
                    tabuList[nextTabu0] = i + (D / 150) + random.randint(1, 10)
                    currentSolution = nextSolution
                    currentFitness = nextFitness
                    fitnessList = fitList

                # print(currentFitness)
                # 更新最佳值
                if currentFitness < self.BestFitness:
                    self.BestSolution = currentSolution
                    self.BestFitness = currentFitness
                    bestList = fitnessList
                    # print('best updated: ')
                    # print(self.BestSolution)
                    # print(self.BestFitness)

            # 停机
            if total_iter > self._maxRunningTicks:
                return

            total_iter += 1

            # 跳坑
            if bestList is not None:
                bestBits = np.argpartition(-1 *
                                           bestList, kth=oneThird)[:oneThird]
                # print('flip:', bestBits)
                # print('from', currentSolution.reshape([1, -1]), 'jump to:')
                currentSolution[bestBits] = 1 - currentSolution[bestBits]
                # print(currentSolution.reshape([1, -1]))
                currentFitness = self._fit(currentSolution)
                # print(str(self.BestFitness) + ' -> ' + str(currentFitness))
                tabuList[:] = 0
                tabuList[bestBits] = (D / 150) + random.randint(1, 10)


class ITS_Instance():
    def __init__(self, maxIteration: int, maxRunningTicks: int, C, query):
        """禁忌搜索初始化

        Args:
            maxIteration (int): tabu search 迭代次数
            maxRunningTicks (int): 停机时间, 秒
        """

        self.BestSolution = None
        self.BestFitness = 0.0

        self._maxIteration = maxIteration
        self._maxRunningTicks = maxRunningTicks

        self._C = C
        self._x = query

        self._k = C.shape[0]

    def _fit(self, b):
        # faster than matrix operation, see test_norm.py
        return np.linalg.norm(b.dot(self._C) - self._x)

    def _matfit(self, b):
        # norm for each row, [N, K] * [K, D] - [N, D]
        return np.linalg.norm(b.dot(self._C) - self._x, axis=1)

    def _initSolution(self):
        b = np.zeros([self._k, ], dtype=np.int8)
        """ Primal test, half of b set to 1, and another half is 0 """
        choice = np.random.choice(self._k, self._k // 2)
        b[choice] = 1
        return b

    def _updateFitnessList(self, solution):
        # get the fitness list that using this solution and flip each bit
        D = solution.shape[0]
        # [D, D]
        neighbors = np.tile(solution, [D, 1])
        # flip bits for every row
        neighbors[range(D), range(D)] = 1 - neighbors[range(D), range(D)]

        # compute fit function for each row
        return self._matfit(neighbors)

    def run(self):

        currentSolution = self._initSolution()

        D = currentSolution.shape[0]
        oneThird = D // 4

        currentFitness = self._fit(currentSolution)

        tabuList = np.zeros([D], dtype=np.float)

        self.BestFitness = currentFitness
        self.BestSolution = currentSolution

        # 非禁忌候选集中的解
        nextSolution = None
        nextFitness = int(sys.maxsize)
        nextTabu0 = None

        # 禁忌候选集中的解
        tabuSolution = None
        netxtTabuFitness = int(sys.maxsize)
        nextTabu1 = None

        total_iter = 0

        while True:
            bestList = None
            tabuList = np.zeros([D], dtype=np.float)
            # 迭代限制
            for i in range(self._maxIteration):
                nextFitness = int(sys.maxsize)
                netxtTabuFitness = int(sys.maxsize)

                success = False

                # 邻域候选集中查找

                # [D, D]
                neighbors = np.tile(currentSolution, [D, 1])
                # flip bits for every row
                neighbors[range(D), range(D)] = 1 - neighbors[range(D), range(D)]

                # compute fit function for each row
                vals = self._matfit(neighbors)

                for j in range(D):
                    # 是否满足禁忌要求
                    if i > tabuList[j]:
                        # 选择邻域中评价最好的解作为候选
                        if nextFitness > vals[j]:
                            nextSolution = neighbors[j]
                            nextFitness = vals[j]
                            # 记录位置
                            nextTabu0 = j
                            success = True
                    else:
                        # 选择邻域中评价最好的解作为候选
                        if netxtTabuFitness > vals[j]:
                            tabuSolution = neighbors[j]
                            netxtTabuFitness = vals[j]
                            # 记录位置
                            nextTabu1 = j

                # 赦免：禁忌表已满，或禁忌中的值比最好值还好
                if not success or (nextFitness > netxtTabuFitness and netxtTabuFitness < self.BestFitness):
                    # 更新禁忌表，当前值
                    # tabuList[nextTabu1] = i + (D / 150) + random.randint(1, 10)
                    tabuList[:] = 0
                    tabuList[nextTabu1] = i + (D / 150) + random.randint(1, 10)
                    currentSolution = tabuSolution
                    currentFitness = netxtTabuFitness
                else:
                    # 更新禁忌表，当前值
                    tabuList[nextTabu0] = i + (D / 150) + random.randint(1, 10)
                    currentSolution = nextSolution
                    currentFitness = nextFitness

                # print(currentFitness)
                # 更新最佳值
                if currentFitness < self.BestFitness:
                    self.BestSolution = currentSolution
                    self.BestFitness = currentFitness
                    bestList = self._updateFitnessList(self.BestSolution)
                    # print('best updated: ')
                    # print(self.BestSolution)
                    # print(self.BestFitness)

            # 停机
            if total_iter > self._maxRunningTicks:
                return

            total_iter += 1

            # 跳坑
            if bestList is not None:
                bestBits = np.argpartition(bestList, kth=oneThird)[:oneThird]
                # print('flip:', bestBits)
                # print('from', currentSolution.reshape([1, -1]), 'jump to:')
                currentSolution[bestBits] = 1 - currentSolution[bestBits]
                # print(currentSolution.reshape([1, -1]))
                currentFitness = self._fit(currentSolution)
                # print(str(self.BestFitness) + ' -> ' + str(currentFitness))
                tabuList[:] = 0
                tabuList[bestBits] = (D / 150) + random.randint(1, 10)
<file_sep>/my_common.py
# encoding: utf-8
import numpy as np
import math
import random

#根据给定的范围随机选取规定个数的数字
def random_select(left, right, candidate):
    return random.sample(range(left, right), candidate)

#根据给定范围随机将部分值的符号改变
def random_change_sgn(left, right, candidate, vec):
    list = random_select(left, right, candidate)
    for i in list:
        vec[i] = -vec[i]
    return vec

#根据大于0则为1，其他为0为0的原则对向量进行二值化
def binazation(vec):
    return np.int64(vec > 0.5)

#随机一个长度为len的01向量
def random_vecb(len):
    return np.array([i % 2 for i in np.random.randint(1, 100, size=(len))])

def random_matB(row,col):
    matB = []
    for i in range(col):
        matB = np.hstack((matB, random_vecb(row)))
    return np.array(matB).reshape(col, row).T

#计算两个矩阵间的距离
def matrix_distance(mat, norm=1):
    if norm == 1:
        fmat = np.fabs(mat)
        return np.sum(fmat) / mat.size
    else:
        fmat = np.multiply(mat, mat)
        return np.sqrt(np.sum(fmat))/mat.size

#输入的x为标准矩阵，处理时需要转置
def trans_intmat_binmat(x, hd):
    ture_x = x.T
    row = ture_x.shape[0]
    col = ture_x.shape[1]
    binmat = []

    for i in range(row):
        binvec = trans_int10_vector2(ture_x[i], hd)
        binmat = np.hstack((binmat, binvec))

    return np.array(binmat).reshape(row, col*hd).T

def trans_b10_vector2(int_b,hd, is_v=True):
    #math.floor(x)向下取整
    #横向并联两个矩阵或向量np.hstack((a,b))
    b_vector = np.zeros(0)
    int_b = int_b.squeeze()
    for num in int_b:
        b_vector = np.hstack((b_vector, trans_int10_vector2(math.floor(num),hd)))
    if is_v:
        return b_vector.T
    else:
        return b_vector

def trans_int10_vector2(x, hd):
    num = x
    bin_vector = np.zeros(hd)
    while True:
        hd = hd - 1
        num, remainder = divmod(num, 2)
        if remainder == 1:
            bin_vector[hd] = 1
        if hd == 0:
            return bin_vector

#输入的x为标准矩阵，处理时需要先转置
def trans_binmat_intmat(x, hd):
    ture_x = x.T
    row = ture_x.shape[0]
    col = ture_x.shape[1]
    intmat = []

    for i in range(row):
        intvec = trans_binvec_intvec(ture_x[i], hd)
        intmat = np.hstack((intmat, intvec))

    return np.array(intmat).reshape(row, col//hd).T

def trans_binvec_intvec(x, hd):
    len = x.shape[0]
    mul = 0
    intvec = np.zeros(len//hd)
    while True:
        vec = x[hd * mul:hd * (mul+1)]
        intvec[mul] = trans_binvec_int(vec)
        mul = mul + 1
        if hd * mul == len:
            return intvec

def trans_binvec_int(x):
    len = x.shape[0]
    num = 0
    for i in range(len):
        num = num + x[i] * math.pow(2, len - i - 1)
    return num

def up2dim(x):
    x = np.expand_dims(x, axis=2)
    x = np.expand_dims(x, axis=3)
    return x









<file_sep>/random_c.py
import numpy as np

class RandomC(object):
    def __init__(self, X, seed=4):
        self.r = np.random
        self.r.seed(seed)
        self.x_min = X.min()
        self.x_max = X.max()
        self.x_mean = X.mean()
        self.x_std = X.std()
        self.x_d = X.shape[0]
        self.c = []

    def product_c(self, k):
        if k % 3 == 0:
            sub_k = k // 3
            self.c = np.hstack((self.uniform_c(sub_k), self.randn_c(sub_k), self.rand_c(sub_k)))
        elif k % 2 == 0:
            sub_k = k // 2
            self.c = np.hstack((self.uniform_c(sub_k), self.randn_c(sub_k)))
        else:
            self.c = self.uniform_c(k)
        return self.c

    def uniform_c(self, sub_k):
        return self.r.uniform(self.x_min, self.x_max, size=(self.x_d, sub_k))

    def randn_c(self, sub_k):
        return self.x_std * self.r.randn(self.x_d, sub_k) + self.x_mean

    def rand_c(self, sub_k):
        return self.r.rand(self.x_d, sub_k) * self.x_mean
