Skip to content

Commit f8ae14a

Browse files
professorcode1professorcode1
committed
metropolis hasting done and broken
1 parent b1cf3fc commit f8ae14a

File tree

11 files changed

+284
-38
lines changed

11 files changed

+284
-38
lines changed

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Application.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ static inline bool file_exists(const std::string& name) {
1111
const unsigned int SCREEN_WIDTH = 1920;
1212
const unsigned int SCREEN_HEIGHT = 1024;
1313
Application::Application() :screen{
14-
[this](boost::gregorian::date date, float value) {
14+
[this](boost::posix_time::ptime date, float value) {
1515
this->data.insert(std::make_pair(date, value));
1616
//std::cout << boost::gregorian::to_simple_string(this->data.begin()->first) << std::endl;;
1717
},
@@ -35,7 +35,7 @@ Application::Application() :screen{
3535
dataIterator++;
3636
nextDataIterator++;
3737
}
38-
this->HMC_Wiggins = new HMC(this->parameter, std::move(StocksData), this->deviceNameToWorkload);
38+
this->HMC_Wiggins = new MetropolisHasting(this->parameter, std::move(StocksData), this->deviceNameToWorkload);
3939
},
4040
[this]() {
4141
this->HMC_Wiggins->iterate();
@@ -72,6 +72,7 @@ Application::Application() :screen{
7272
},
7373
SCREEN_WIDTH, SCREEN_HEIGHT
7474
} {
75+
std::cout << "Hello??" << std::endl;
7576
glfwSetErrorCallback(glfw_error_callback);
7677
if (!glfwInit()){
7778
std::cout << "glfwInit failed. Terminating program" << std::endl;

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Application.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,15 @@
1212
#include "vendor/stbi/stbi.h"
1313
#include "Engine/SYCLComputer.h"
1414
#include "Engine/HMC.h"
15+
#include "Engine/MetropolisHasting.h"
1516

1617

1718
class Application
1819
{
1920
private:
2021
GLFWwindow* window;
2122
Screen screen;
22-
std::map<boost::gregorian::date, float> data;
23+
std::map<boost::posix_time::ptime, float> data;
2324
std::map<std::string, int> deviceNameToWorkload;
2425
AlgorithmParameter parameter;
2526
WigginsAlgorithm* HMC_Wiggins;

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/Algorithm.cpp‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ WigginsAlgorithm::WigginsAlgorithm(
1313
parameter.m_volatility_sigma.generateBuffer(parameter.m_DiscretCountOfContinuiosSpace)
1414
}
1515
{
16+
DumpToFileFloatVector(stocksdata);
1617
this->iteration_count = 0;
1718
int workloadTotal = 0;
1819
std::vector<cl::sycl::device> devices = cl::sycl::device::get_devices();
@@ -172,3 +173,26 @@ AlgorithmResponse WigginsAlgorithm::get_response() {
172173
return this->m_response;
173174
}
174175

176+
177+
178+
extern SYCL_EXTERNAL float probability(
179+
const ProbabilityDomain domain,
180+
int timeSeriesLength,
181+
oneapi::dpl::minstd_rand& engine,
182+
const cl::sycl::accessor<float, 1, sycl::access::mode::read>& returns
183+
) {
184+
const float minus_log_p = potential_energy_U(domain, timeSeriesLength, engine, returns);
185+
return exp(-1.0f * minus_log_p);
186+
}
187+
188+
uint32_t DistributionIndex(float lower, float upper, uint32_t divisions, float value) {
189+
float delx = (upper - lower) / divisions;
190+
uint32_t res = floor((value - lower) / delx);
191+
if (res >= divisions) {
192+
return divisions - 1;
193+
}
194+
if (res < 0) {
195+
return 0;
196+
}
197+
return res;
198+
}

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/Algorithm.h‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,7 @@ class WigginsAlgorithm
2929
ProbabilityDomain start;
3030

3131
virtual void BurnIn() = 0;
32-
/*cl::sycl::event run_inference(
33-
cl::sycl::queue& q, uint32_t number_of_simulations,
34-
uint32_t number_of_days, cl::sycl::buffer<float, 2> data);*/
32+
3533
public:
3634
WigginsAlgorithm(
3735
const AlgorithmParameter &parameter,
@@ -73,4 +71,13 @@ extern SYCL_EXTERNAL ProbabilityDomain potential_energy_gradient(
7371
const AlgorithmParameter& parameters,
7472
oneapi::dpl::minstd_rand& engine,
7573
const cl::sycl::accessor<float, 1, sycl::access::mode::read>& returns
76-
);
74+
);
75+
76+
extern SYCL_EXTERNAL float probability(
77+
const ProbabilityDomain domain,
78+
int timeSeriesLength,
79+
oneapi::dpl::minstd_rand& engine,
80+
const cl::sycl::accessor<float, 1, sycl::access::mode::read>& returns
81+
);
82+
83+
uint32_t DistributionIndex(float lower, float upper, uint32_t divisions, float value);

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/HMC.cpp‎

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -104,17 +104,8 @@ void HMC::BurnInInternal() {
104104
this->start.sigma = resultHostAccessor[2];
105105
}
106106

107-
uint32_t DistributionIndex(float lower, float upper, uint32_t divisions, float value) {
108-
float delx = (upper - lower) / divisions;
109-
uint32_t res = floor((value - lower) / delx);
110-
if (res == divisions) {
111-
return res - 1;
112-
}
113-
return res;
114-
}
115-
116107

117-
cl::sycl::event exectue_wiggins_algorithm_on_device(
108+
static cl::sycl::event exectue_wiggins_algorithm_on_device(
118109
cl::sycl::queue& q, const AlgorithmDeviceData& device_data, int workItems,
119110
const AlgorithmParameter& parameter, const ProbabilityDomain& start,
120111
cl::sycl::buffer<float, 1> returns, cl::sycl::buffer<float, 1>& theta,

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/MetropolisHasting.cpp‎

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,201 @@
11
#include "MetropolisHasting.h"
2+
3+
MetropolisHasting::MetropolisHasting(
4+
const AlgorithmParameter& parameter,
5+
const std::vector<float>& stocksdata,
6+
const std::map<std::string, int>& deviceNameToWorkload
7+
) :WigginsAlgorithm{ parameter ,stocksdata , deviceNameToWorkload } {
8+
this->BurnIn();
9+
}
10+
11+
static int sycl_strogest_cpu_selector(const sycl::device& d) {
12+
return d.is_cpu();
13+
}
14+
15+
bool mh_sample(
16+
const ProbabilityDomain current_q,
17+
ProbabilityDomain& new_domain,
18+
const AlgorithmParameter& parameters,
19+
oneapi::dpl::normal_distribution<float>& theta_normal,
20+
oneapi::dpl::normal_distribution<float>& mu_normal,
21+
oneapi::dpl::normal_distribution<float>& sigma_normal,
22+
oneapi::dpl::uniform_real_distribution<float>& zero_one_uniform,
23+
int tsLength,
24+
oneapi::dpl::minstd_rand& engine,
25+
const cl::sycl::accessor<float, 1, sycl::access::mode::read>& returns
26+
) {
27+
ProbabilityDomain p{ theta_normal(engine) , mu_normal(engine), sigma_normal(engine) };
28+
ProbabilityDomain q = current_q + p;
29+
const float current_f = probability(current_q, tsLength, engine, returns);
30+
const float purposed_f = probability(q, tsLength, engine, returns);
31+
const bool accepted = (purposed_f / current_f) > zero_one_uniform(engine);
32+
new_domain = ((int)accepted) * q + ((int)(1 ^ accepted)) * current_q;
33+
return accepted;
34+
}
35+
36+
void MetropolisHasting::BurnInInternal() {
37+
using namespace cl::sycl;
38+
queue q(sycl_strogest_cpu_selector);
39+
unsigned int seed = std::chrono::system_clock::now().time_since_epoch().count();
40+
buffer<float, 1> result(range<1>{3});
41+
q.submit([&, this](handler& h) {
42+
const AlgorithmParameter parameter = this->m_parameter;
43+
accessor returnsAccessor(m_returns, h, read_only);
44+
accessor resultAccessor(result, h, write_only);
45+
int tsLength = m_returns.size();
46+
ProbabilityDomain start = this->start;
47+
h.single_task([=]() {
48+
oneapi::dpl::minstd_rand engine(seed, 0);
49+
oneapi::dpl::normal_distribution<float> theta_normal(0.0f, parameter.m_volatility_theta.guassian_step_sd);
50+
oneapi::dpl::normal_distribution<float> mu_normal(0.0f, parameter.m_volatility_mu.guassian_step_sd);
51+
oneapi::dpl::normal_distribution<float> sigma_normal(0.0f, parameter.m_volatility_sigma.guassian_step_sd);
52+
oneapi::dpl::uniform_real_distribution<float> zero_one_uniform(0.0f, 1.0f);
53+
ProbabilityDomain oldDomain{ start }, nextDomain;
54+
for (int i = 0; i < parameter.m_BurnIn; i++) {
55+
mh_sample(
56+
oldDomain,
57+
nextDomain,
58+
parameter,
59+
theta_normal,
60+
mu_normal,
61+
sigma_normal,
62+
zero_one_uniform,
63+
tsLength,
64+
engine,
65+
returnsAccessor
66+
);
67+
oldDomain = nextDomain;
68+
}
69+
resultAccessor[0] = nextDomain.theta;
70+
resultAccessor[1] = nextDomain.mu;
71+
resultAccessor[2] = nextDomain.sigma;
72+
});
73+
}).wait();
74+
host_accessor resultHostAccessor(result, read_write);
75+
this->start.theta = resultHostAccessor[0];
76+
this->start.mu = resultHostAccessor[1];
77+
this->start.sigma = resultHostAccessor[2];
78+
}
79+
80+
static cl::sycl::event exectue_wiggins_algorithm_on_device(
81+
cl::sycl::queue& q, const AlgorithmDeviceData& device_data, int workItems,
82+
const AlgorithmParameter& parameter, const ProbabilityDomain& start,
83+
cl::sycl::buffer<float, 1> returns, cl::sycl::buffer<float, 1>& theta,
84+
cl::sycl::buffer<float, 1>& mu, cl::sycl::buffer<float, 1>& sigma
85+
) {
86+
using namespace cl::sycl;
87+
88+
int tsLength = returns.size();
89+
unsigned int seed = std::chrono::system_clock::now().time_since_epoch().count();
90+
//run the MNC Algorithm, write result to theta mu sigma buffer
91+
return q.submit([&](handler& h) {
92+
accessor thetaAccessor(theta, h, write_only);
93+
accessor muAccessor(mu, h, write_only);
94+
accessor sigmaAccessor(sigma, h, write_only);
95+
accessor returnsAccessor(returns, h, read_only);
96+
h.parallel_for(range<1>(workItems), [=](id<1> thread_id) {
97+
oneapi::dpl::minstd_rand engine(seed, thread_id);
98+
oneapi::dpl::normal_distribution<float> theta_normal(0.0f, parameter.m_volatility_theta.guassian_step_sd);
99+
oneapi::dpl::normal_distribution<float> mu_normal(0.0f, parameter.m_volatility_mu.guassian_step_sd);
100+
oneapi::dpl::normal_distribution<float> sigma_normal(0.0f, parameter.m_volatility_sigma.guassian_step_sd);
101+
oneapi::dpl::uniform_real_distribution<float> zero_one_uniform(0.0f, 1.0f);
102+
ProbabilityDomain oldDomain{ start }, nextDomain;
103+
mh_sample(
104+
oldDomain,
105+
nextDomain,
106+
parameter,
107+
theta_normal,
108+
mu_normal,
109+
sigma_normal,
110+
zero_one_uniform,
111+
tsLength,
112+
engine,
113+
returnsAccessor
114+
);
115+
thetaAccessor[thread_id] = nextDomain.theta;
116+
muAccessor[thread_id] = nextDomain.mu;
117+
sigmaAccessor[thread_id] = nextDomain.sigma;
118+
});
119+
});
120+
121+
}
122+
123+
void MetropolisHasting::iterateInternal() {
124+
/*
125+
using namespace cl::sycl;
126+
std::vector<cl::sycl::event> events;
127+
std::vector<cl::sycl::buffer<float, 1> > thetaHists;
128+
std::vector<cl::sycl::buffer<float, 1> > muHists;
129+
std::vector<cl::sycl::buffer<float, 1> > sigmaHists;
130+
float mu_lower = this->m_parameter.m_volatility_mu.mean - this->m_parameter.m_volatility_mu.buffer_range_sigma_multiplier * this->m_parameter.m_volatility_mu.sd;
131+
float mu_higher = this->m_parameter.m_volatility_mu.mean + this->m_parameter.m_volatility_mu.buffer_range_sigma_multiplier * this->m_parameter.m_volatility_mu.sd;
132+
uint32_t DiscreteQuanta = this->m_parameter.m_DiscretCountOfContinuiosSpace;
133+
for (int deviceIndex = 0; deviceIndex < this->m_QueuesAndData.size(); deviceIndex++) {
134+
const int workItems = ceil(this->m_parameter.m_GraphUpdateIteration * this->m_QueuesAndData[deviceIndex].second.m_workload_fraction);
135+
136+
thetaHists.emplace_back(workItems);
137+
muHists.emplace_back(workItems);
138+
sigmaHists.emplace_back(workItems);
139+
events.push_back(
140+
exectue_wiggins_algorithm_on_device(
141+
this->m_QueuesAndData[deviceIndex].first,
142+
this->m_QueuesAndData[deviceIndex].second, workItems,
143+
this->m_parameter, this->start,
144+
this->m_returns,
145+
thetaHists.back(), muHists.back(), sigmaHists.back()
146+
)
147+
);
148+
}
149+
for (int deviceIndex = 0; deviceIndex < this->m_QueuesAndData.size(); deviceIndex++) {
150+
const int workItems = ceil(this->m_parameter.m_GraphUpdateIteration * this->m_QueuesAndData[deviceIndex].second.m_workload_fraction);
151+
host_accessor thetaList(thetaHists[deviceIndex], read_only);
152+
for (int index = 0; index < workItems; index++) {
153+
float theta = thetaList[index];
154+
uint32_t thetaHistIndex = DistributionIndex(
155+
this->m_parameter.m_volatility_theta.lower,
156+
this->m_parameter.m_volatility_theta.upper,
157+
DiscreteQuanta,
158+
theta
159+
);
160+
this->m_response.theta.data[thetaHistIndex] += 1;
161+
this->m_response.theta.sum += theta;
162+
}
163+
164+
host_accessor muList(muHists[deviceIndex], read_only);
165+
for (int index = 0; index < workItems; index++) {
166+
float mu = muList[index];
167+
uint32_t muHistIndex = DistributionIndex(
168+
mu_lower,
169+
mu_higher,
170+
DiscreteQuanta,
171+
mu
172+
);
173+
this->m_response.mu.data[muHistIndex] += 1;
174+
this->m_response.mu.sum += mu;
175+
176+
}
177+
178+
179+
host_accessor sigmaList(sigmaHists[deviceIndex], read_only);
180+
for (int index = 0; index < workItems; index++) {
181+
float sigma = sigmaList[index];
182+
uint32_t thetaHistIndex = DistributionIndex(
183+
this->m_parameter.m_volatility_sigma.lower,
184+
this->m_parameter.m_volatility_sigma.upper,
185+
this->m_parameter.m_DiscretCountOfContinuiosSpace,
186+
sigma
187+
);
188+
this->m_response.theta.data[thetaHistIndex] += 1;
189+
this->m_response.theta.sum += sigma;
190+
}
191+
}
192+
*/
193+
this->iteration_count++;
194+
}
195+
196+
void MetropolisHasting::BurnIn() {
197+
this->processing = std::async(std::launch::async, &MetropolisHasting::BurnInInternal, this);
198+
}
199+
void MetropolisHasting::iterate() {
200+
this->processing = std::async(std::launch::async, &MetropolisHasting::iterateInternal, this);
201+
}

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/MetropolisHasting.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,14 +4,15 @@ class MetropolisHasting:public WigginsAlgorithm
44
{
55
private:
66
void BurnIn() override;
7+
void BurnInInternal();
78
public:
89
MetropolisHasting(
910
const AlgorithmParameter& parameter,
1011
const std::vector<float>& stocksdata,
1112
const std::map<std::string, int>& deviceNameToWorkload
1213
);
13-
14-
14+
15+
void iterateInternal();
1516
void iterate() override;
1617

1718
};

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/Utilities.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,4 +52,20 @@ float evaluate_normal(float x, float mean, float sd) {
5252
x /= sd;
5353
x /= sqrtf(2 * M_PI);
5454
return x;
55+
}
56+
57+
58+
void DumpToFileFloatVector(
59+
std::vector<float> data
60+
) {
61+
std::ofstream file("returns dump.csv");
62+
if (file.is_open()) {
63+
for (const auto d : data) {
64+
file << d<<std::endl;
65+
}
66+
file.close();
67+
}
68+
else {
69+
throw std::runtime_error("Cannot open file to dump returns");
70+
}
5571
}

‎src/MonteCarloStockPrediction/MonteCarloStockPrediction/Engine/Utilities.h‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,9 @@ std::vector<float> ReadCSVExtractFloatColumns(
1414

1515
float geometric_mean(float a, float r, int n);
1616

17-
float evaluate_normal(float x, float mean, float sd);
17+
float evaluate_normal(float x, float mean, float sd);
18+
19+
20+
void DumpToFileFloatVector(
21+
std::vector<float> data
22+
);

0 commit comments

Comments
 (0)