Skip to content

Instantly share code, notes, and snippets.

@hxhc

hxhc/CUSim.cpp

Last active May 6, 2022 11:31
Show Gist options
  • Save hxhc/2cb8fc583ea137b7b2653048792168f3 to your computer and use it in GitHub Desktop.
Save hxhc/2cb8fc583ea137b7b2653048792168f3 to your computer and use it in GitHub Desktop.
Download ZIP
CUSim
Raw
bind.cpp
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <pybind11/stl.h>
#include "CUSim.hpp"
// 编译器链接时需要提供python3.lib
// 运行时需要将python3.dll放入exe文件根目录中
namespace py = pybind11;
PYBIND11_MODULE(cusim, m) {
m.def("sim_sd", &sim_sd, "simulate sd");
m.def("sim_mu", &sim_mu, "simulate mu");
}
Raw
CUSim.cpp
#include "CUTest.hpp"
#include <omp.h>
std::vector<Eigen::ArrayXd> sim_sd(double mu, double start, double end, int sd_num=61, int sim_num = 10000, double low_lvl = 85.0, double upp_lvl = 115.0,
double ck1 = 2.2, double ck2 = 0.25, double ck3 = 1.7, double cL = 15.0, double uk1 = 2.4, double uk2 = 2.0, double uL1 = 15.0, double uL2 = 25.0) {
int batch = 20000; // simulate 20k tablets in one batch
Eigen::ArrayXd failrate = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_fail_s1 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_fail_s2 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_pass_s1 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_pass_s2 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd usp_fail_s1 = Eigen::ArrayXd::Zero(sd_num); // no fail in stage 1 test indeed
Eigen::ArrayXd usp_fail_s2 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd usp_pass_s1 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd usp_pass_s2 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_acc_s1 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd chp_acc = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd usp_acc_s1 = Eigen::ArrayXd::Zero(sd_num);
Eigen::ArrayXd usp_acc = Eigen::ArrayXd::Zero(sd_num);
for (int i = 0; i < sd_num; i++) {
double sd = start + (end - start) * i / double(sd_num - 1);
// normal distribution
std::random_device device;
Eigen::Rand::Vmt19937_64 urng{ device() };
Eigen::Rand::NormalGen<double> norm_gen{mu, sd};
Eigen::ArrayXd X = norm_gen.generate<Eigen::ArrayXd>(batch, 1, urng);
failrate[i] = double((X > upp_lvl || X < low_lvl).count()) / batch;
//result = np.array([pass_s1, pass_s2, fail_s1, fail_s2])
Eigen::Array4d chp_result = Eigen::Array4d::Zero(4);
Eigen::Array4d usp_result = Eigen::Array4d::Zero(4);
// random draw 30 samples, first 10 as s1, last 20 as s2
Eigen::Rand::UniformIntGen<int> uni_int_gen{ 0, 19999 };
Eigen::Rand::Vmt19937_64 urng_uni_int{ device() };
Eigen::ArrayXi rng_idx = uni_int_gen.generate<Eigen::ArrayXi>(30, 1, urng_uni_int);
for (int j = 0; j < sim_num; ++j) {
rng_idx = uni_int_gen.generate<Eigen::ArrayXi>(30, 1, urng_uni_int);
Eigen::ArrayXi idx_s1 = rng_idx(Eigen::seq(0, 9));
Eigen::ArrayXi idx_s2 = rng_idx(Eigen::seq(10, 29));
Eigen::Array<double, 10, 1> s1 = X(idx_s1);
Eigen::Array<double, 20, 1> s2 = X(idx_s2);
ChPTest chp = ChPTest(s1, s2, ck1, ck2, ck3, cL);
USPTest usp = USPTest(s1, s2, uk1, uk2, uL1, uL2);
chp_result += chp.get_result();
usp_result += usp.get_result();
}
chp_pass_s1[i] = chp_result[0] / double(sim_num);
chp_pass_s2[i] = chp_result[1] / double(sim_num);
chp_fail_s1[i] = chp_result[2] / double(sim_num);
chp_fail_s2[i] = chp_result[3] / double(sim_num);
usp_pass_s1[i] = usp_result[0] / double(sim_num);
usp_pass_s2[i] = usp_result[1] / double(sim_num);
usp_fail_s1[i] = usp_result[2] / double(sim_num);
usp_fail_s2[i] = usp_result[3] / double(sim_num);
if (failrate[i] <= 0.05) {
chp_acc_s1[i] = chp_pass_s1[i];
chp_acc[i] = chp_pass_s1[i] + chp_pass_s2[i];
usp_acc_s1[i] = usp_pass_s1[i];
usp_acc[i] = usp_pass_s1[i] + usp_pass_s2[i];
}
else {
chp_acc_s1[i] = chp_fail_s1[i];
chp_acc[i] = chp_fail_s1[i] + chp_fail_s2[i];
usp_acc[i] = usp_fail_s2[i];
}
}
std::vector<Eigen::ArrayXd> result{ failrate, chp_acc_s1, chp_acc, usp_acc_s1, usp_acc, chp_pass_s1, chp_pass_s2, usp_pass_s1, usp_pass_s2, usp_fail_s1, usp_fail_s2};
return result;
}
std::vector<Eigen::ArrayXd> sim_mu(double sd, double start, double end, int mu_num, int sim_num, double low_lvl, double upp_lvl,
double ck1, double ck2, double ck3, double cL, double uk1, double uk2, double uL1, double uL2) {
int batch = 20000; // simulate 20k tablets in one batch
Eigen::ArrayXd failrate = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_fail_s1 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_fail_s2 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_pass_s1 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_pass_s2 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd usp_fail_s1 = Eigen::ArrayXd::Zero(mu_num); // no fail in stage 1 test indeed
Eigen::ArrayXd usp_fail_s2 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd usp_pass_s1 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd usp_pass_s2 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_acc_s1 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd chp_acc = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd usp_acc_s1 = Eigen::ArrayXd::Zero(mu_num);
Eigen::ArrayXd usp_acc = Eigen::ArrayXd::Zero(mu_num);
for (int i = 0; i < mu_num; i++) {
double mu = start + (end - start) * i / double(mu_num - 1);
// normal distribution
std::random_device device;
Eigen::Rand::Vmt19937_64 urng{ device() };
Eigen::Rand::NormalGen<double> norm_gen{ mu, sd };
Eigen::ArrayXd X = norm_gen.generate<Eigen::ArrayXd>(batch, 1, urng);
failrate[i] = double((X > upp_lvl || X < low_lvl).count()) / batch;
//result = np.array([pass_s1, pass_s2, fail_s1, fail_s2])
Eigen::Array4d chp_result = Eigen::Array4d::Zero(4);
Eigen::Array4d usp_result = Eigen::Array4d::Zero(4);
// random draw 30 samples, first 10 as s1, last 20 as s2
Eigen::Rand::UniformIntGen<int> uni_int_gen{ 0, 19999 };
Eigen::Rand::Vmt19937_64 urng_uni_int{ device() };
Eigen::ArrayXi rng_idx = uni_int_gen.generate<Eigen::ArrayXi>(30, 1, urng_uni_int);
for (int j = 0; j < mu_num; ++j) {
rng_idx = uni_int_gen.generate<Eigen::ArrayXi>(30, 1, urng_uni_int);
Eigen::ArrayXi idx_s1 = rng_idx(Eigen::seq(0, 9));
Eigen::ArrayXi idx_s2 = rng_idx(Eigen::seq(10, 29));
Eigen::Array<double, 10, 1> s1 = X(idx_s1);
Eigen::Array<double, 20, 1> s2 = X(idx_s2);
ChPTest chp = ChPTest(s1, s2, ck1, ck2, ck3, cL);
USPTest usp = USPTest(s1, s2, uk1, uk2, uL1, uL2);
chp_result += chp.get_result();
usp_result += usp.get_result();
}
chp_pass_s1[i] = chp_result[0] / double(mu_num);
chp_pass_s2[i] = chp_result[1] / double(mu_num);
chp_fail_s1[i] = chp_result[2] / double(mu_num);
chp_fail_s2[i] = chp_result[3] / double(mu_num);
usp_pass_s1[i] = usp_result[0] / double(mu_num);
usp_pass_s2[i] = usp_result[1] / double(mu_num);
usp_fail_s1[i] = usp_result[2] / double(mu_num);
usp_fail_s2[i] = usp_result[3] / double(mu_num);
if (failrate[i] <= 0.05) {
chp_acc_s1[i] = chp_pass_s1[i];
chp_acc[i] = chp_pass_s1[i] + chp_pass_s2[i];
usp_acc_s1[i] = usp_pass_s1[i];
usp_acc[i] = usp_pass_s1[i] + usp_pass_s2[i];
}
else {
chp_acc_s1[i] = chp_fail_s1[i];
chp_acc[i] = chp_fail_s1[i] + chp_fail_s2[i];
usp_acc[i] = usp_fail_s2[i];
}
}
std::vector<Eigen::ArrayXd> result{ failrate, chp_acc_s1, chp_acc, usp_acc_s1, usp_acc, chp_pass_s1, chp_pass_s2, usp_pass_s1, usp_pass_s2, usp_fail_s1, usp_fail_s2 };
return result;
}
Raw
CUSim.hpp
#pragma once
#include "CUTest.hpp"
/// @brief simulate CUtest under different standard deviation values and a fixed mean value
/// @param mu mean value of the normal distribution
/// @param start the start sd
/// @param end the end sd
/// @param mu_num the nunmber of sd values to be generated, defalt is 101
/// @param sim_num the simulation number, default is 10000
/// @param low_lvl the lower level of qualified tablets
/// @param upp_lvl the upper level of qualified tablets
/// @param ck1 the value of k1 in ChPtest
/// @param ck2 the value of k2 in ChPtest
/// @param ck3 the value of k3 in ChPtest
/// @param cL the value of L in ChPtest
/// @param uk1 the value of k1 in USPtest
/// @param uk2 the value of k2 in USPtest
/// @param uL1 the value of L1 in USPtest
/// @param uL2 the value of L2 in USPtest
/// @return vector of [failrate, chp_acc_s1, chp_acc, usp_acc_s1, usp_acc, chp_pass_s1, chp_pass_s2, usp_pass_s1, usp_pass_s2, usp_fail_s1, usp_fail_s2]
std::vector<Eigen::ArrayXd> sim_sd(double mu, double start, double end, int sd_num, int sim_num, double low_lvl, double upp_lvl,
double ck1, double ck2, double ck3, double cL, double uk1, double uk2, double uL1, double uL2);
/// @brief simulate CUtest under different mean values and a fixed standard deviation
/// @param mu standard deviation of the normal distribution
/// @param start the start mu
/// @param end the end mu
/// @param mu_num the nunmber of mu values to be generated, defalt is 101
/// @param sim_num the simulation number, default is 10000
/// @param low_lvl the lower level of qualified tablets
/// @param upp_lvl the upper level of qualified tablets
/// @param ck1 the value of k1 in ChPtest
/// @param ck2 the value of k2 in ChPtest
/// @param ck3 the value of k3 in ChPtest
/// @param cL the value of L in ChPtest
/// @param uk1 the value of k1 in USPtest
/// @param uk2 the value of k2 in USPtest
/// @param uL1 the value of L1 in USPtest
/// @param uL2 the value of L2 in USPtest
/// @return vector of [failrate, chp_acc_s1, chp_acc, usp_acc_s1, usp_acc, chp_pass_s1, chp_pass_s2, usp_pass_s1, usp_pass_s2, usp_fail_s1, usp_fail_s2]
std::vector<Eigen::ArrayXd> sim_mu(double sd, double start, double end, int mu_num, int sim_num, double low_lvl, double upp_lvl,
double ck1, double ck2, double ck3, double cL, double uk1, double uk2, double uL1, double uL2);
Raw
CUTest.cpp
#include "CUTest.hpp"
// standard deviation with 1 degree of freedom (unbiased)
double std_ddof(Eigen::ArrayXd array, int ddof=1) {
double standard_deviation = std::sqrt((array - array.mean()).square().sum() / (array.size() - ddof));
return standard_deviation;
}
ChPTest::ChPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in) {
// check length error
if (s1_in.size() != 10) {
throw("The length of s1 must be 10!");
}
if (s2_in.size() != 20) {
throw("The legnth of s2 must be 20!");
}
s1 = s1_in;
s2 = s2_in;
stage_2 << s1, s2;
mean_s1 = s1.mean();
std_s1 = std_ddof(s1, 1);
mean_stage2 = stage_2.mean();
std_stage2 = std_ddof(stage_2, 1);
}
ChPTest::ChPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in, const double k1_in, const double k2_in, const double k3_in, const double L_in) {
// check length error
if (s1_in.size() != 10) {
throw("The length of s1 must be 10!");
}
if (s2_in.size() != 20) {
throw("The legnth of s2 must be 20!");
}
s1 = s1_in;
s2 = s2_in;
k1 = k1_in;
k2 = k2_in;
k3 = k3_in;
L = L_in;
stage_2 << s1, s2;
mean_s1 = s1.mean();
std_s1 = std_ddof(s1, 1);
mean_stage2 = stage_2.mean();
std_stage2 = std_ddof(stage_2, 1);
}
ChPTest::~ChPTest() { }
//defintion of member function get_result
Eigen::Array4d ChPTest::get_result() {
int pass_s1 = 0;
int pass_s2 = 0;
int fail_s1 = 0;
int fail_s2 = 0;
double Ave = std::abs(100 - mean_s1);
if (Ave + k1 * std_s1 <= L) {
pass_s1 = 1;
}
else if (Ave + std_s1 > L) {
fail_s1 = 1;
}
else {
Ave = std::abs(100 - mean_stage2);
if (Ave <= k2 * L) {
if (pow(Ave, 2) + pow(std_stage2, 2) <= k2 * pow(L, 2)) {
pass_s2 = 1;
}
else {
fail_s2 = 1;
}
}
else
{
if (Ave + k3 * std_stage2 <= L) {
pass_s2 = 1;
}
else {
fail_s2 = 1;
}
}
}
Eigen::Array4d result;
result << pass_s1, pass_s2, fail_s1, fail_s2;
return result;
}
USPTest::USPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in) {
// check length error
if (s1_in.size() != 10) {
throw("The length of s1 must be 10!");
}
if (s2_in.size() != 20) {
throw("The legnth of s2 must be 20!");
}
s1 = s1_in;
s2 = s2_in;
stage_2 << s1, s2;
mean_s1 = s1.mean();
std_s1 = std_ddof(s1, 1);
mean_stage2 = stage_2.mean();
std_stage2 = std_ddof(stage_2, 1);
}
USPTest::USPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in, double k1_in, double k2_in, double L1_in, double L2_in)
{
// check length error
if (s1_in.size() != 10) {
throw("The length of s1 must be 10!");
}
if (s2_in.size() != 20) {
throw("The legnth of s2 must be 20!");
}
s1 = s1_in;
s2 = s2_in;
k1 = k1_in;
k2 = k2_in;
L1 = L1_in;
L2 = L2_in;
stage_2 << s1, s2;
mean_s1 = s1.mean();
std_s1 = std_ddof(s1, 1);
mean_stage2 = stage_2.mean();
std_stage2 = std_ddof(stage_2, 1);
}
USPTest::~USPTest() {}
// return Eigen::Array4d [pass_s1, pass_s2, fail_s1, fail_s2]
Eigen::Array4d USPTest::get_result() {
int pass_s1 = 0;
int pass_s2 = 0;
int fail_s1 = 0; //no fail in stage 1, indeed
int fail_s2 = 0;
double Ave = 0;
double Ave2 = 0;
double limit = 0;
if (mean_s1 < 98.5) {
Ave = 98.5 - mean_s1;
}
else if (mean_s1 > 101.5) {
Ave = mean_s1 - 101.5;
}
else {
Ave = 0;
}
// stage 1 test
if (Ave + k1 * std_s1 <= L1) {
pass_s1 = 1;
}
else {
if (mean_stage2 < 98.5) {
Ave2 = 98.5 - mean_stage2;
limit = 98.5;
}
else if (mean_stage2 > 101.5) {
Ave2 = mean_stage2 - 101.5;
limit = 101.5;
}
else {
Ave2 = 0;
limit = mean_stage2;
}
if ((Ave2 + k2 * std_stage2 <= L1) && (stage_2 >= (1 - 0.01 * L2) * limit).all() && (stage_2 <= (1 + 0.01 * L2) * limit).all()) {
pass_s2 = 1;
}
else {
fail_s2 = 1;
}
}
Eigen::Array4d result;
result << pass_s1, pass_s2, fail_s1, fail_s2;
return result;
}
Raw
CUtest.hpp
#pragma once
#include <Eigen/Dense>
#include <EigenRand/EigenRand>
#include <chrono>
#include <cmath>
#include <random>
#include <iostream>
/// @brief standard deviation with 1 degree of freedom (unbiased)
/// @param array an Eigen array
/// @param ddof the degree of freedom, default is 1
/// @return the standard deviation of the array
double std_ddof(Eigen::ArrayXd array, int ddof);
class ChPTest {
public:
// members
Eigen::Array<double, 10, 1> s1;
Eigen::Array<double, 20, 1> s2;
Eigen::Array<double, 30, 1> stage_2;
double k1{ 2.2 };
double k2{ 0.5 };
double k3{ 1.7 };
double L{ 15 };
double std_s1;
double std_stage2;
double mean_s1;
double mean_stage2;
Eigen::Array4d get_result();
// constructor function
ChPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in);
ChPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in, const double k1_in, const double k2_in, const double k3_in, const double L_in);
~ChPTest();
};
class USPTest {
public:
Eigen::Array<double, 10, 1> s1;
Eigen::Array<double, 20, 1> s2;
Eigen::Array<double, 30, 1> stage_2;
double k1{ 2.4 };
double k2{ 2.0 };
double L1{ 15.0 };
double L2{ 25.0 };
double std_s1;
double std_stage2;
double mean_s1;
double mean_stage2;
// constructor function
USPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in);
USPTest(const Eigen::ArrayXd& s1_in, const Eigen::ArrayXd& s2_in, double k1_in, double k2_in, double L1_in, double L2_in);
~USPTest();
// member function
Eigen::Array4d get_result();
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment