# 数模模板

### 层次分析法

import numpy as np
import pandas as pd
import warnings
class AHP:
def __init__(self, criteria, a):
self.RI = (0, 0, 0.528, 0.93, 1.973, 1.24, 1.532, 1.41, 21.45, 12.49)
self.criteria = criteria
self.b = b
self.num_criteria = criteria.shape[0]
self.num_project = b[0].shape[0]
def cal_weights(self, input_matrix):
input_matrix = np.array(input_matrix)
n, n1 = input_matrix.shape
assert n == n1, 'ERROR'
for i in range(n):
for j in range(n):
if np.abs(input_matrix[i, j] * input_matrix[j, i] - 1) > 1e-7:
raise ValueError('ERROR')
eigenvalues, eigenvectors = np.linalg.eig(input_matrix)
max_idx = np.argmax(eigenvalues)
max_eigen = eigenvalues[max_idx].real
eigen = eigenvectors[:, max_idx].real
eigen = eigen / eigen.sum()
if n > 9:
CR = None
warnings.warn('ERROR')
else:
CI = (max_eigen - n) / (n - 1)
CR = CI / self.RI[n]
print(CR)
return max_eigen, CR, eigen
def run(self):
max_eigen, CR, criteria_eigen = self.cal_weights(self.criteria)
print('MAX{:<5f},CR={:<5f},EXAM{}PASS'.format(max_eigen, CR, '' if CR < 0.1 else 'no'))
print('WEIGHT={}\n'.format(criteria_eigen))

max_eigen_list, CR_list, eigen_list = [], [], []
for i in self.b:
max_eigen, CR, eigen = self.cal_weights(i)
max_eigen_list.append(max_eigen)
CR_list.append(CR)
eigen_list.append(eigen)

pd_print = pd.DataFrame(eigen_list,
index=['ZZ' + str(i) for i in range(self.num_criteria)],
columns=['PLAN' + str(i) for i in range(self.num_project)],
)
pd_print.loc[:, 'MAX'] = max_eigen_list
pd_print.loc[:, 'CR'] = CR_list
pd_print.loc[:, 'TOTAL'] = pd_print.loc[:, 'CR'] < 0.1
print(pd_print)
obj = np.dot(criteria_eigen.reshape(1, -1), np.array(eigen_list))
print('The best plan is :{}'.format(np.argmax(obj)))
return obj

if __name__ == '__main__':
criteria = np.array([[1,2,3,5,4],
[1/2,1,5,2,4],
[1/3,1/5,1,1/3,1/4],
[1/5,1/2,3,1,1/2],
[1/4,1/4,4,2,1]])
b1 = np.array([[1, 1 / 3, 1 / 8], [3, 1, 1 / 3], [8, 3, 1], [8, 3, 1], [8, 3, 1]])
b2 = np.array([[1, 2, 5], [1 / 2, 1, 2], [1 / 5, 1 / 2, 1], [8, 3, 1], [8, 3, 1]])
b = [b1, b2]
a = AHP(criteria, b).run()


### TOPsis

import math
import csv
f_best=np.array([0.40848,0.22385,0.48667,0.42196])
f_worst=np.array([0.01225,0.00284,0.00002,0.00036])
f = open(r"E:\mathorcup\处理结果.csv",'w',newline="")
csv_writer = csv.writer(f)
csv_writer.writerow(['材料名称', '密度', '弹性模量', '产量', '拉伸强度','负理想解', '正理想解','最终得分'])
list=data.values.tolist()
list.remove(['材料名称', '密度', '弹性模量', '产量', '拉伸强度'])
for i in list:
f=np.array([float(i[1]),float (i[2]),float(i[3]),float(i[4])])
csv_writer.writerow(i)


### 热卡

import pandas as pd
import math
import operator
import csv
from math import radians, cos, sin, asin, sqrt
list_point=[]
dict={}
list=data.values.tolist()
f = open(r"e:\停车点替换点.csv",'w',newline="")
csv_writer = csv.writer(f)
csv_writer.writerow(["停车点坐标","出发需求","到达需求","需求总数","可替代停车点坐标"])
def geodistance(lng1,lat1,lng2,lat2):
lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)])
dlon=lng2-lng1
dlat=lat2-lat1
a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
distance=2*asin(sqrt(a))*6371*1000
distance=round(distance,3)
return distance
for i in range(len(list)):
point_temp=list[i][0].replace("(","").replace(")","").split(",")
if point_temp not in list_point:
list_point.append(point_temp)

for i in range(len(list)):
list_write=list[i]
if(int(list[i][3])<10):
point_temp=list[i][0].replace("(","").replace(")","").split(",")
flag=0
for j in list_point:
if(flag==0):
distance=geodistance(point_temp[0],point_temp[1],j[0],j[1])
if(distance<150 and distance!=0):
list_write.append("("+str(j[0])+","+str(j[1])+")")
csv_writer.writerow(list_write)
flag=1
continue
if(flag==0):
list_write.append("无替代点")
csv_writer.writerow(list_write)
else:
list_write.append("无需去除")
csv_writer.writerow(list_write)
f.close()



### 小波模型

# -*- coding: utf-8 -*-
import numpy as np
import math
import matplotlib.pyplot as plt
import pandas as pd
import datetime
from scipy import interpolate
from pandas import DataFrame,Series

import numpy as np
import pywt

data = np.linspace(1, 4, 7)

# pywt.threshold方法讲解：
#               pywt.threshold（data，value，mode ='soft'，substitute = 0 ）
#               data：数据集，value：阈值，mode：比较模式默认soft，substitute：替代值，默认0，float类型

#data:   [ 1.   1.5  2.   2.5  3.   3.5  4. ]
#output：[ 6.   6.   0.   0.5  1.   1.5  2. ]
#soft 因为data中1小于2，所以使用6替换，因为data中第二个1.5小于2也被替换，2不小于2所以使用当前值减去2，，2.5大于2，所以2.5-2=0.5.....

print(pywt.threshold(data, 2, 'soft',6))

#data:   [ 1.   1.5  2.   2.5  3.   3.5  4. ]
#hard data中绝对值小于阈值2的替换为6，大于2的不替换
print (pywt.threshold(data, 2, 'hard',6))

#data:   [ 1.   1.5  2.   2.5  3.   3.5  4. ]
#data中数值小于阈值的替换为6，大于等于的不替换
print (pywt.threshold(data, 2, 'greater',6) )

print (data  )
#data:   [ 1.   1.5  2.   2.5  3.   3.5  4. ]
#data中数值大于阈值的，替换为6
print (pywt.threshold(data, 2, 'less',6) )

# 算法竞赛

### 判断小于MAXN的数是不是素数

/*
*  素数筛选，判断小于MAXN的数是不是素数
*  notprime是一张表，false表示是素数，true表示不是
*/

const int MAXN = 1000010;
bool notprime[MAXN];

void init()
{
memset(notprime, false, sizeof(notprime));
notprime[0] = notprime[1] = true;
for (int i = 2; i < MAXN; i++)
{
if (!notprime[i])
{
if (i > MAXN / i)   //  阻止后边i * i溢出（或者i,j用long long)
{
continue;
}
//  直接从i * i开始就可以，小于i倍的已经筛选过了
for (int j = i * i; j < MAXN; j += i)
{
notprime[j] = true;
}
}
}
}

### 生成连续素数表

/*
*  素数筛选，查找出小于等于MAXN的素数
*  prime[0]存素数的个数
*/

const int MAXN = 100000;
int prime[MAXN + 1];

void getPrime()
{
memset(prime, 0, sizeof(prime));
for (int i = 2; i <= MAXN; i++)
{
if (!prime[i])
{
prime[++prime[0]] = i;
}
for (int j = 1; j <= prime[0] && prime[j] <= MAXN / i; j++)
{
prime[prime[j] * i] = 1;
if (i % prime[j] == 0)
{
break;
}
}
}
}

### 组合数C(n, r)

int com(int n, int r)   //  return C(n, r)
{
if (n - r > r)
{
r = n - r;      //  C(n, r) = C(n, n - r)
}
int i, j, s = 1;
for (i = 0, j = 1; i < r; ++i)
{
s *= (n - i);
for (; j <= r && s % j == 0; ++j)
{
s /= j;
}
}
return s;
}

### 集合划分问题

/*
*  n元集合分划为k类的方案数记为S(n, k),称为第二类Stirling数。
*  如{A,B,C}可以划分{{A}, {B}, {C}}, {{A, B}, {C}}, {{B, C}, {A}}, {{A, C}, {B}}, {{A, B, C}}。
*  即一个集合可以划分为不同集合(1...n个)的种类数
*  CALL: compute(N); 每当输入一个n,输出B[n]
*/
const int N = 2001;
int data[N][N], B[N];

void NGetM(int m, int n)    // m 个数 n 个集合
{
//  data[i][j]: i个数分成j个集合
int min, i, j;
data[0][0] = 1;
for (i = 1; i <= m; i++)
{
data[i][0] = 0;
}
for (i = 0; i <= m; i++)
{
data[i][i + 1] = 0;
}
for (i = 1; i <= m; i++)
{
if (i < n)
{
min = i;
}
else
{
min = n;
}
for (j = 1; j <= min; j++)
{
data[i][j] = (j * data[i - 1][j] + data[i - 1][j - 1]);
}
}
return ;
}

void compute(int m)
{
//  b[i]: Bell数
NGetM(m, m);
memset(B, 0, sizeof(B));
int i, j;
for (i=1; i <= m; i++)
{
for (j = 0; j <= i; j++)
{
B[i] += data[i][j];
}
}
return ;
}

### 基姆拉尔森公式

/*
*  已知1752年9月3日是Sunday，并且日期控制在1700年2月28日后
*/
char name[][15] = { "monday", "tuesday", "wednesday", "thursday", "friday", "saturday", "sunday"};

int main()
{
int d, m, y, a;
printf("Day: ");
scanf("%d", &d);
printf("Month: ");
scanf("%d", &m);
printf("Year: ");
scanf("%d", &y);
//  1月2月当作前一年的13,14月
if (m == 1 || m == 2)
{
m += 12;
y--;
}
//  判断是否在1752年9月3日之前,实际上合并在一起倒更加省事
if ((y < 1752) || (y == 1752 && m < 9) || (y == 1752 && m == 9 && d < 3))
{
//  因为日期控制在1700年2月28日后，所以不用考虑整百年是否是闰年
a = (d + 2 * m + 3 * (m + 1) / 5 + y + y / 4 + 5) % 7;
}
else
{
//  这里需要考虑整百年是否是闰年的情况
a = (d + 2 * m + 3 * (m + 1) / 5 + y + y / 4 - y / 100 + y / 400) % 7;  //  实际上这个可以当做公式背下来
}
printf("it's a %s\n", name[a]);
return 0;
}

### 矩阵原理单独求解（斐波那契）

/*
*  求斐波那契数列第N项，模MOD
*/
#define mod(a, m) ((a) % (m) + (m)) % (m)

const int MOD = 1e9 + 9;

struct MATRIX
{
long long a[2][2];
};

MATRIX a;
long long f[2];

void ANS_Cf(MATRIX a)
{
f[0] = mod(a.a[0][0] + a.a[1][0], MOD);
f[1] = mod(a.a[0][1] + a.a[1][1], MOD);
return ;
}

MATRIX MATRIX_Cf(MATRIX a, MATRIX b)
{
MATRIX ans;
int k;
for (int i = 0; i < 2; i++)
{
for (int j = 0; j < 2; j++)
{
ans.a[i][j] = 0;
k = 0;
while (k < 2)
{
ans.a[i][j] += a.a[k][i] * b.a[j][k];
ans.a[i][j] = mod(ans.a[i][j], MOD);
++k;
}
}
}
return ans;
}

MATRIX MATRIX_Pow(MATRIX a, long long n)
{
MATRIX ans;
ans.a[0][0] = 1;
ans.a[1][1] = 1;
ans.a[0][1] = 0;
ans.a[1][0] = 0;
while (n)
{
if (n & 1)
{
ans = MATRIX_Cf(ans, a);
}
n = n >> 1;
a = MATRIX_Cf(a, a);
}
return ans;
}

int main()
{
long long n;
while (cin >> n)
{
if (n == 1)
{
cout << '1' << '\n';
continue;
}
a.a[0][0] = a.a[0][1] = a.a[1][0] = 1;
a.a[1][1] = 0;
a = MATRIX_Pow(a, n - 2);
ANS_Cf(a);
cout << f[0] << '\n';
}
return 0;
}

### AC自动机

/*
*  求目标串中出现了几个模式串
*/
struct Trie
{
int next[500010][26], fail[500010], end[500010];
int root, L;
int newnode()
{
for (int i = 0; i < 26; i++)
{
next[L][i] = -1;
}
end[L++] = 0;
return L - 1;
}

void init()
{
L = 0;
root = newnode();
}

void insert(char buf[])
{
int len = (int)strlen(buf);
int now = root;
for (int i = 0; i < len; i++)
{
if (next[now][buf[i] - 'a'] == -1)
{
next[now][buf[i] - 'a'] = newnode();
}
now = next[now][buf[i] - 'a'];
}
end[now]++;
}

void build()
{
queue<int>Q;
fail[root] = root;
for (int i = 0; i < 26; i++)
{
if (next[root][i] == -1)
{
next[root][i] = root;
}
else
{
fail[next[root][i]] = root;
Q.push(next[root][i]);
}
}
while (!Q.empty())
{
int now = Q.front();
Q.pop();
for (int i = 0;i < 26;i++)
{
if (next[now][i] == -1)
{
next[now][i] = next[fail[now]][i];
}
else
{
fail[next[now][i]]=next[fail[now]][i];
Q.push(next[now][i]);
}
}
}
}

int query(char buf[])
{
int len = (int)strlen(buf);
int now = root;
int res = 0;
for (int i = 0; i < len; i++)
{
now = next[now][buf[i] - 'a'];
int temp = now;
while (temp != root)
{
res += end[temp];
end[temp] = 0;
temp = fail[temp];
}
}
return res;
}

void debug()
{
for (int i = 0; i < L; i++)
{
printf("id = %3d,fail = %3d,end = %3d,chi = [", i, fail[i], end[i]);
for (int j = 0; j < 26; j++)
{
printf("%2d", next[i][j]);
}
printf("]\n");
}
}
};

char buf[1000010];
Trie ac;

int main()
{
int T;
int n;
scanf("%d", &T);
while(T--)
{
scanf("%d", &n);
ac.init();
for (int i = 0; i < n; i++)
{
scanf("%s", buf);
ac.insert(buf);
}
ac.build();
scanf("%s", buf);
printf("%d\n", ac.query(buf));
}
return 0;
}

### Prim

/*
*  Prim求MST
*  耗费矩阵cost[][]，初始化为INF，标号从0开始，0 ~ n－1
*  返回最小生成树的权值，返回-1表示原图不连通
*/

const int INF = 0x3f3f3f3f;
const int MAXN = 110;
bool vis[MAXN];
int lowc[MAXN];
int cost[MAXN][MAXN];

//  修正cost（添加边）
void updata(int x, int y, int v)
{
cost[x - 1][y - 1] = v;
cost[y - 1][x - 1] = v;
return ;
}

int Prim(int cost[][MAXN], int n)   //  0 ~ n - 1
{
int ans = 0;
memset(vis, false, sizeof(vis));
vis[0] = true;
for (int i = 1; i < n; i++)
{
lowc[i] = cost[0][i];
}
for (int i = 1; i < n; i++)
{
int minc = INF;
int p = -1;
for (int j = 0; j < n; j++)
{
if (!vis[j] && minc > lowc[j])
{
minc = lowc[j];
p = j;
}
}
if (minc == INF)
{
return -1;  //  原图不连通
}
ans += minc;
vis[p] = true;
for (int j = 0; j < n; j++)
{
if (!vis[j] && lowc[j] > cost[p][j])
{
lowc[j] = cost[p][j];
}
}
}
return ans;
}

### Kruskal

/*
*  Kruskal算法求MST
*  对边操作，并排序
*  切记：初始化赋值问题（tol）
*/

const int MAXN = 110;   //  最大点数
const int MAXM = 10000; //  最大边数

int F[MAXN];    //  并查集使用

struct Edge
{
int u;      //  起点
int v;      //  终点
int w;      //  权值
} edge[MAXM];   //  存储边的信息

int tol;        //  边数，加边前赋值为0

void addEdge(int u, int v, int w)
{
edge[tol].u = u;
edge[tol].v = v;
edge[tol++].w = w;
return ;
}

bool cmp(Edge a, Edge b)
{
//  排序函数，将边按照权值从小到大排序
return a.w < b.w;
}

int find(int x)
{
if (F[x] == x)
{
return x;
}
else
{
return F[x] = find(F[x]);
}
}

int Kruskal(int n)  //  传入点数，返回最小生成树的权值，如果不连通则返回-1
{
for (int i = 0; i <= n; i++)
{
F[i] = i;
}
sort(edge, edge + tol, cmp);

int cnt = 0;    //  计算加入的边数
int ans = 0;
for (int i = 0; i < tol; i++)
{
int u = edge[i].u;
int v = edge[i].v;
int w = edge[i].w;
int tOne = find(u);
int tTwo = find(v);
if (tOne != tTwo)
{
ans += w;
F[tOne] = tTwo;
cnt++;
}
if (cnt == n - 1)
{
break;
}
}
if (cnt < n - 1)
{
return -1;  //  不连通
}
else
{
return ans;
}
}

### 并查集

/*
*  INIT: makeset(n);
*  CALL: findset(x); unin(x, y);
*/
const int N = 1010;

struct lset
{
int p[N], rank[N], sz;
{
if (x == y)
{
return ;
}
if (rank[x] > rank[y])
{
p[y] = x;
}
else
{
p[x] = y;
}
if (rank[x] == rank[y])
{
rank[y]++;
}
return ;
}
void makeset(int n)
{
sz = n;
for (int i = 0; i < sz; i++)
{
p[i] = i;
rank[i] = 0;
}
return ;
}
int findset(int x)
{
if (x != p[x])
{
p[x] = findset(p[x]);
}
return p[x];
}
void unin(int x, int y)
{
};