mathorcup2023年第三问 创新解法完整程序

完整程序

import pandas as pd
import re
import numpy as np
import pandas as pd
from dwave.system import DWaveSampler, LeapHybridSampler, EmbeddingComposite
from pyqubo import Array, Binary
import json
import multiprocessing as mp
import sys

def partition(i):
    # Read credit card data
    data = pd.read_csv("附件1：data_100.csv", header=0)
    t = data.iloc[:, ::2] # Approval rate
    h = data.iloc[:, 1::2]  # Bad debt rate
    # Calculate matrix shape
    _, n_cols = t.shape
    # Calculate column ranges for each partition
    t1_cols = [i]  # Initialize as a list with a single element
    if i+1 <= n_cols-1 and i+50 <= n_cols-1:
        t2_cols = list(range(i+1, i+51))
    elif i+1 <= n_cols-1 and i+50 > n_cols-1:
        t2_cols = list(range(i+1, n_cols)) + list(range(0, 50-(n_cols-i)))
    elif i == n_cols-1:
        t2_cols = list(range(0, 50))
    t3_cols = list(set(range(n_cols)) - set(t1_cols) - set(t2_cols))
    t1 = t.iloc[:, t1_cols].values
    t2 = t.iloc[:, t2_cols].values
    t3 = t.iloc[:, t3_cols].values
     # Partition h
    h1 = h.iloc[:, t1_cols].values
    h2 = h.iloc[:, t2_cols].values
    h3 = h.iloc[:, t3_cols].values
     # Print results
    print("t1 shape: ", t1.shape)
    print("t2 shape: ", t2.shape)
    print("t3 shape: ", t3.shape)
    print("h1 shape: ", h1.shape)
    print("h2 shape: ", h2.shape)
    print("h3 shape: ", h3.shape)
    print(t1.T.reshape(1, -1).shape)
    print(t2.T.reshape(1, -1).shape)
    print(t3.T.reshape(1, -1).shape)
    return t1.T.reshape(1, -1) ,t2.T.reshape(1, -1),t3.T.reshape(1, -1),h1.T.reshape(1, -1),h2.T.reshape(1, -1),h3.T.reshape(1, -1)

def get_result(solution,j):
    # x, y 用于存储下标
    x1, y1 = [], []
    x2, y2 = [], []
    x3, y3 = [], []
    i = 0
    pattern_idx = re.compile(r"[0-9]+")
    for key, value in solution.items():
        i += 1
        print(i)
        if i >1000:
             break
        
        if i <=10 and value == 1:
                print(key, value)
                _,x_idx, y_idx = pattern_idx.findall(key)
                print(f"第{j+1}张卡选择第{int(y_idx)+1}个阈值")
                x1.append(j)
                y1.append(int(y_idx))
               
        elif i > 10 and i<=510  and value == 1:
                print(key, value)
                _,_, y_idx = pattern_idx.findall(key)
                #int(y_idx)+1是因为y_idx从0开始，而阈值从1开始，(int(y_idx)+1)//10+1是哪怕是等于0，也要加1
                a=(int(y_idx)+1)//10+1
                b=int(y_idx)%10
                if a+j+1<=100:
                    print(f"第{a+j+1}张卡选择第{b+1}个阈值")
                    x2.append(a+j)
                    y2.append(b)
                else:
                    print(f"第{a+j+1-100}张卡选择第{b+1}个阈值")
                    x2.append(a+j-100)
                    y2.append(b)                
                
        elif i >= 510 and i<1000 and value == 1:
                print(key, value)
                _,_, y_idx = pattern_idx.findall(key)
                a=(int(y_idx)+1)//10+1
                b=int(y_idx)%10
                if j+1+50+a<=100:
                    print(f"第{j+1+50+a}张卡选择第{b+1}个阈值")
                    x3.append(j+50+a)
                    y3.append(b)
                else:
                    print(f"第{j+1+50+a-100}张卡选择第{b+1}个阈值")
                    x3.append(j+50+a-100)
                    y3.append(b)                
    print(x1[0],y1[0])
    print(x2[0],y2[0])
    print(x3[0],y3[0])
    C = 1000000  # 总贷款资金
    R = 0.08     # 利率
    n = 100      # 信用评分卡数量
    data = pd.read_csv("附件1：data_100.csv", header=0)
    t = data.iloc[:, ::2] # Approval rate
    h = data.iloc[:, 1::2]  # Bad debt rate           
    # T, H 分别表示总通过率和总坏账率
    T = t.iloc[y1[0], x1[0]] * t.iloc[y2[0], x2[0]] * t.iloc[y3[0], x3[0]]
    H = 1 / 3 * (h.iloc[y1[0], x1[0]] + h.iloc[y2[0], x2[0]] + h.iloc[y3[0], x3[0]])
    result = C * T * (R * (1 - H) - H)
    print(f"result\n")
    return result,x1,y1,x2,y2,x3,y3

def operation(i):
    # 定义常数和变量
    C = 1000000  # 总贷款资金
    R = 0.08     # 利率
    n = 100      # 信用评分卡数量
    P = C        # 惩罚系数

    # 定义目标函数
    y1 = Array.create('y1', shape=(1, 10), vartype='BINARY')
    y2 = Array.create('y2', shape=(1, 500), vartype='BINARY')
    y3 = Array.create('y3', shape=(1, 490), vartype='BINARY')
    z = Array.create('z', shape=(10, 490), vartype='BINARY')

    t1,t2,t3,h1,h2,h3=partition(i)

    Q = 0

    for i in range(10):
        for j in range(500):
            for k in range(490):
                T_ijk = t1[0, i] * t2[0, j] * t3[0, k]
                H_ijk = (h1[0, i] + h2[0, j] + h3[0, k]) / 3
                Q += y2[0, j] * z[i, k] * T_ijk * (C * R * (1 - H_ijk) - C * H_ijk)

    c1 = 99999999 * (sum(y1[0,:]) - 1)**2 
    
    c2 = 10000000 * (sum(y2[0,:]) - 1)**2 
    
    c3 = 15000000 * (sum(y3[0,:]) - 1)**2 

    c4 = 0

    for i in range(10):
        for k in range(490):
            c4 += 100000000 * (3 * z[i, k] - 2 * z[i, k] * y1[0, i] -
                    2 * z[i, k] * y3[0, k] + y1[0, i] * y3[0, k])

    model = -Q + c1 + c2 + c3 + c4
    model = model.compile()
    Q, _ = model.to_qubo()

    # 求解 QUBO 模型
    sampler = LeapHybridSampler()        
    # 使用 LeapHybridSampler
    response = sampler.sample_qubo(Q)   
    print("最优解:")
    solution = response.first.sample
    return solution


def process_operation(i):
    a = []
    solution = operation(i)
    result, x1, y1, x2, y2, x3, y3 = get_result(solution, i)
    result_dict = {'result': result, 'x1': x1, 'y1': y1,
                   'x2': x2, 'y2': y2, 'x3': x3, 'y3': y3}
    a.append({i: result_dict})
    with open('result.json', 'a') as f:
        json.dump(a, f, ensure_ascii=False)
        f.write('\n')

if __name__ == '__main__':
    sys.stdout = open('/dev/null', 'w')
    pool = mp.Pool(processes=64) # Use 4 processes for concurrent execution
    for i in range(0,99):
        pool.apply_async(process_operation, args=(i,))
    pool.close()
    pool.join()

    max_result = float('-inf')  # Initialize max_result to negative infinity
    max_row = None  # Initialize max_row to None
    with open('result.json', 'r') as f:
        for line in f:
            row = json.loads(line)[0]  # Extract the dictionary from the list
            for key, value in row.items():
                result = value['result']
                if result > max_result:
                    max_result = result
                    max_row = row
    # Output the row with the highest 'result' value
    print(f"最大值是: {max_row[next(iter(max_row))]['result']}")
    print(f"第一张信用卡取第: {max_row[next(iter(max_row))]['x1'][0]+1}张, 选择 {max_row[next(iter(max_row))]['y1'][0]+1}作为阈值:")
    print(f"第二张信用卡取第: {max_row[next(iter(max_row))]['x2'][0]+1}张, 选择 {max_row[next(iter(max_row))]['y2'][0]+1}作为阈值:")
    print(f"第三张信用卡取第: {max_row[next(iter(max_row))]['x3'][0]+1}张, 选择 {max_row[next(iter(max_row))]['y3'][0]+1}作为阈值:")