Skip to content

Instantly share code, notes, and snippets.

@2bbb

2bbb/logistic_regression.hpp

Last active January 8, 2016 09:40
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save 2bbb/571ab176dbe7cb500bba to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save 2bbb/571ab176dbe7cb500bba to your computer and use it in GitHub Desktop.
Download ZIP
Raw
logistic_regression.hpp
//
// logistic_regression.hpp
//
// Created by ISHII 2bit on 2016/01/08.
//
//
#pragma once
#include "vec.hpp"
#include <vector>
namespace logistic_regression {
using default_value_type = float;
template <std::size_t num, typename value_type = default_value_type>
struct model {
using record = bbb::vec<num, value_type>;
using result = bool;
struct datum {
record xs;
result y;
};
using data = std::vector<datum>;
using offset = value_type;
using weight = bbb::vec<num, value_type>;
struct parameter {
offset w0;
weight ws;
};
static inline std::int8_t value(const datum &m) {
return m.y ? 1 : -1;
}
static inline value_type sigmoid(value_type t) {
return 1.0 / (1.0 + exp(-t));
}
static inline value_type prediction(const parameter &p, const datum &m) {
return p.w0 + p.ws * m.xs;
}
static inline value_type probability(const parameter &p, const datum &m) {
return sigmoid(prediction(p, m));
}
static inline value_type error(const parameter &p, const datum &m) {
return log(1 + exp(-value(m) * prediction(p, m)));
}
static inline value_type objective(const parameter &p, const data &d, value_type smooth_factor = 0.01) {
value_type sum{(p.ws * p.ws) * smooth_factor};
std::for_each(d.begin(), d.end(), [&](const datum &m){ sum += error(p, m); });
return sum;
}
static inline value_type d_w0(const parameter &p, const datum &m) {
value_type tmp = -value(m) * prediction(p, m);
return -value(m) * exp(tmp) / (1 + exp(tmp));
}
static inline value_type d_ws(const parameter &p, const datum &m, std::size_t index, value_type smooth_factor = 0.01) {
value_type tmp = -value(m) * prediction(p, m);
return m.xs[index] * -value(m) * exp(tmp) / (1 + exp(tmp)) + 2 * p.ws[index] * smooth_factor;
}
static std::vector<std::size_t> random_indicies;
static inline parameter step(const parameter &p, const data &d, float scale = 0.01f) {
if(random_indicies.empty()) {
for(std::size_t i = 0; i < d.size(); i++) random_indicies.push_back(i);
std::random_shuffle(random_indicies.begin(), random_indicies.end());
}
if(d.size() <= random_indicies.back()) {
random_indicies.pop_back();
return p;
}
const datum &m = d[random_indicies.back()];
random_indicies.pop_back();
parameter p0{p};
p0.w0 -= scale * d_w0(p, m);
for(std::size_t i = 0; i < p.ws.size(); i++) p0.ws[i] -= scale * d_ws(p, m, i);
return p0;
}
};
template <std::size_t num, typename value_type>
std::vector<std::size_t> model<num, value_type>::random_indicies;
};
Raw
ofApp.cpp
#include "ofMain.h"
#include "logistic_regression.hpp"
class ofApp : public ofBaseApp {
using model = logistic_regression::model<2>;
model::data data;
model::parameter p;
ofEasyCam cam;
float factor;
float correctness;
std::vector<model::datum> targets;
public:
void setup() {
data.clear();
targets.clear();
factor = 0.025;
correctness = -1;
p = {ofRandom(1.0f), {ofRandom(1.0f), ofRandom(1.0f)}};
ofLogNotice("param::weights") << p.ws[0] << ", " << p.ws[1];
ofVec2f p0{ofRandom(0.0f, 1.0f), ofRandom(0.0f, 1.0f)}, q0{ofRandom(0.0f, 1.0f), ofRandom(0.0f, 1.0f)};
for(std::size_t i = 0; i < 1000; i++) {
model::datum m;
m.y = ofRandom(1.0) < 0.6;
float theta = ofRandom(2 * M_PI);
ofVec2f dp = ofRandom(0.2f) * ofVec2f(cos(theta), sin(theta));
m.xs[0] = (m.y ? p0.x : q0.x) + dp.x;
m.xs[1] = (m.y ? p0.y : q0.y) + dp.y;
if(ofRandom(1.0) < 0.1) m.y ^= true;
data.push_back(m);
}
}
void update() {
p = model::step(p, data, factor);
// if(new_correctness / correctness < 0.9) {
// factor = ofClamp(factor * 0.9, 0.0, 0.05);
// correctness = new_correctness;
// } else if(1.1 < new_correctness / correctness){
// factor = ofClamp(factor * 1.1, 0.0, 0.05);
// correctness = new_correctness;
// }
ofLogNotice("correctness") << correctness;
if(0.5 < correctness) {
float new_correctness = model::objective(p, data);
factor = ofClamp(factor * new_correctness / correctness, 0.0, 0.05);
correctness = new_correctness;
} else if(correctness < 0.0) {
correctness = model::objective(p, data);;
} else {
factor = 0.0f;
}
}
void draw() {
ofBackground(0, 0, 0);
for(auto &m : data) {
ofSetColor(model::probability(p, m) < 0.5 ? ofColor::blue : ofColor::green);
if(m.y && model::probability(p, m) < 0.5) ofSetColor(ofColor::red);
if(!m.y && 0.5 <= model::probability(p, m)) ofSetColor(ofColor::red);
if(ofVec2f(m.xs[0], m.xs[1]).distance(ofVec2f(ofGetMouseX() / 720.0f, ofGetMouseY() / 720.0f)) < 0.01f) {
ofDrawBitmapString(ofToString(model::probability(p, m)), 720 * m.xs[0], 720 * m.xs[1] - 16);
}
if(m.y) {
ofDrawCircle(720 * m.xs[0], 720 * m.xs[1], 4);
} else {
ofDrawRectangle(720 * m.xs[0] - 4, 720 * m.xs[1] - 4, 8, 8);
}
}
ofSetColor(255, 0, 0);
ofVec2f sp{0.0f, 0.0f}, ep{1.0f, 1.0f};
if(fabs(p.ws[0]) < fabs(p.ws[1])) {
sp.y = -p.w0 / p.ws[1];
ep.y = -(p.w0 + p.ws[0]) / p.ws[1];
} else {
sp.x = -p.w0 / p.ws[0];
ep.x = -(p.w0 + p.ws[1]) / p.ws[0];
}
sp *= 720.0f;
ep *= 720.0f;
ofDrawLine(sp.x, sp.y, ep.x, ep.y);
for(const auto &m : targets) {
ofSetColor(model::probability(p, m) < 0.5 ? ofColor::blue : ofColor::green);
ofDrawLine((m.xs[0] - 0.01) * 720, m.xs[1] * 720,
(m.xs[0] + 0.01) * 720, m.xs[1] * 720);
ofDrawLine(m.xs[0] * 720, (m.xs[1] - 0.01) * 720,
m.xs[0] * 720, (m.xs[1] + 0.01) * 720);
}
}
void exit() {}
void keyPressed(int key) {
if(key == 'r') {
setup();
}
}
void keyReleased(int key) {}
void mouseMoved(int x, int y) {}
void mouseDragged(int x, int y, int button) {}
void mousePressed(int x, int y, int button) {
model::datum m;
m.xs[0] = x / 720.0f;
m.xs[1] = y / 720.0f;
targets.push_back(m);
}
void mouseReleased(int x, int y, int button) {}
void mouseEntered(int x, int y) {}
void mouseExited(int x, int y) {}
void windowResized(int w, int h) {}
void dragEvent(ofDragInfo dragInfo) {}
void gotMessage(ofMessage msg) {}
};
int main() {
ofSetupOpenGL(720, 720, OF_WINDOW);
ofRunApp(new ofApp);
}
Raw
vec.hpp
//
// vec.hpp
//
// Created by ISHII 2bit on 2016/01/08.
//
//
#pragma once
#include <iostream>
#include <type_traits>
#include <array>
#include <cmath>
#if 201402L <= __cplusplus
# define constexpr_14 constexpr
#else
# define constexpr_14
#endif
namespace bbb {
template <std::size_t s, typename value_t = double>
struct base_vec {
static_assert(std::is_arithmetic<value_t>::value, "required: value_t is arithmetic type");
using value_type = value_t;
base_vec()
: data() {}
base_vec(const base_vec &v)
: data(v.data) {}
template <typename argument, typename ... arguments, typename std::enable_if<std::is_arithmetic<argument>::value>::type *_ = nullptr>
base_vec(argument arg, arguments ... args)
: data({static_cast<value_t>(arg), static_cast<value_t>(args) ...}) {}
template <std::size_t s_, typename value_t_>
base_vec(const base_vec<s_, value_t_> &v) {
for(std::size_t i = 0, end = std::min(s, s_); i < end; i++) data[i] = v[i];
}
constexpr std::size_t size() const { return s; }
value_t &operator[](std::size_t index) { return data[index]; }
constexpr const value_t &operator[](std::size_t index) const { return data[index]; }
value_t &at(std::size_t index) { return data.at(index); }
constexpr const value_t &at(std::size_t index) const { return data.at(index); }
base_vec &operator=(const base_vec &v) { data = v.data; }
base_vec &operator=(const std::array<value_t, s> &v) { data = v; }
template <std::size_t s_, typename value_t_>
base_vec &operator=(const base_vec<s_, value_t_> &v) {
for(std::size_t i = 0, end = std::min(s, s_); i < end; i++) data[i] = v.data[i];
return *this;
}
constexpr const base_vec &operator+() const { return *this; }
constexpr_14 base_vec operator-() const {
base_vec v0;
for(std::size_t i = 0; i < size(); i++) v0 = -data[i];
return v0;
}
constexpr_14 base_vec operator+(const base_vec &v) const {
base_vec v0;
for(std::size_t i = 0; i < size(); i++) v0 = data[i] + v[i];
return v0;
}
base_vec &operator+=(const base_vec &v) {
for(std::size_t i = 0; i < size(); i++) data[i] += v[i];
return *this;
}
constexpr_14 base_vec operator-(const base_vec &v) const {
base_vec v0;
for(std::size_t i = 0; i < size(); i++) v0 = data[i] - v[i];
return v0;
}
base_vec &operator-=(const base_vec &v) {
for(std::size_t i = 0; i < size(); i++) data[i] -= v[i];
return *this;
}
constexpr_14 value_t operator*(const base_vec &v) const {
value_t sum{0.0};
for(std::size_t i = 0; i < size(); i++) sum += data[i] * v[i];
return sum;
}
constexpr value_t dot(const base_vec &v) const {
return *this * v;
}
constexpr_14 base_vec operator*(value_t scale) const {
base_vec v0;
for(std::size_t i = 0; i < size(); i++) v0[i] = data[i] * scale;
return v0;
}
base_vec &operator*=(value_t scale) {
for(std::size_t i = 0; i < size(); i++) data[i] *= scale;
return *this;
}
operator std::array<value_t, s> &() { return data; }
constexpr operator const std::array<value_t, s> &() const { return data; }
constexpr bool operator==(const base_vec &v) const { return data == v.data; }
constexpr bool operator!=(const base_vec &v) const { return data != v.data; }
constexpr bool operator<(const base_vec &v) const { return data < v.data; }
constexpr bool operator<=(const base_vec &v) const { return data <= v.data; }
constexpr bool operator>(const base_vec &v) const { return data > v.data; }
constexpr bool operator>=(const base_vec &v) const { return data >= v.data; }
void swap(base_vec &v) { std::swap(data, v.data); }
using iterator = typename std::array<value_t, s>::iterator;
using const_iterator = typename std::array<value_t, s>::const_iterator;
using reverse_iterator = typename std::array<value_t, s>::reverse_iterator;
using const_reverse_iterator = typename std::array<value_t, s>::const_reverse_iterator;
iterator begin() { return data.begin(); }
iterator end() { return data.end(); }
const_iterator begin() const { return data.cbegin(); }
const_iterator end() const { return data.cend(); }
const_iterator cbegin() const { return data.cbegin(); }
const_iterator cend() const { return data.cend(); }
reverse_iterator rbegin() { return data.rbegin(); }
reverse_iterator rend() { return data.rend(); }
const_reverse_iterator rbegin() const { return data.crbegin(); }
const_reverse_iterator rend() const { return data.crend(); }
const_reverse_iterator crbegin() const { return data.crbegin(); }
const_reverse_iterator crend() const { return data.crend(); }
friend base_vec operator*(value_t scale, const base_vec &v);
protected:
std::array<value_t, s> data;
};
template <std::size_t s, typename value_t = double>
base_vec<s, value_t> operator*(value_t scale, const base_vec<s, value_t> &v) {
return v * s;
}
template <std::size_t s, typename value_t = double>
struct vec : base_vec<s, value_t> {
vec() : base_vec<s, value_t>() {};
template <typename ... arguments>
vec(arguments ... args) : base_vec<s, value_t>(args ...) {};
};
template <typename value_t>
struct vec<2, value_t> : base_vec<2, value_t> {
value_t &x, &y;
vec(value_t x = 0, value_t y = 0)
: base_vec<2, value_t>(x, y)
, x(data[0])
, y(data[1]) {};
template <std::size_t s_, typename value_t_>
vec(const vec<s_, value_t_> &v)
: base_vec<2, value_t>(static_cast<base_vec<s_, value_t>>(v))
, x(data[0])
, y(data[1]) {};
inline vec &operator=(const vec &v) { data = v.data; };
private:
using base_vec<2, value_t>::data;
};
template <typename value_t>
struct vec<3, value_t> : base_vec<3, value_t> {
value_t &x, &y, &z;
vec(value_t x = 0, value_t y = 0, value_t z = 0)
: base_vec<3, value_t>(x, y, z)
, x(data[0])
, y(data[1])
, z(data[2]) {}
operator base_vec<3, value_t> &() { return *this; }
operator const base_vec<3, value_t> &() const { return *this; }
template <std::size_t s_, typename value_t_>
vec(const vec<s_, value_t_> &v)
: base_vec<3, value_t>(static_cast<base_vec<s_, value_t>>(v))
, x(data[0])
, y(data[1])
, z(data[2]) {};
inline vec &operator=(const vec &v) { data = v.data; };
private:
using base_vec<3, value_t>::data;
};
template <typename value_t>
struct vec<4, value_t> : base_vec<4, value_t> {
value_t &x, &y, &z, &w;
vec(value_t x = 0, value_t y = 0, value_t z = 0, value_t w = 0)
: base_vec<4, value_t>(x, y, z, w)
, x(data[0])
, y(data[1])
, z(data[2])
, w(data[3]) {}
template <std::size_t s_, typename value_t_>
vec(const vec<s_, value_t_> &v)
: base_vec<4, value_t>(static_cast<base_vec<s_, value_t>>(v))
, x(data[0])
, y(data[1])
, z(data[2])
, w(data[3]) {};
inline vec &operator=(const vec &v) { data = v.data; };
private:
using base_vec<4, value_t>::data;
};
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment