0418

install_diff/Makefile

+CXX = g++
+BIN = bin
+LIB = lib
+LIB_NAME = autodiff
+OBJS_GRADIENT_DESCENT = root/obj/node.o root/obj/graph.o examples/obj/gradient_descent.o
+OBJS_ANN = root/obj/node.o root/obj/graph.o examples/obj/ann.o
+OBJS_GRADIENT = root/obj/node.o root/obj/graph.o examples/obj/gradient.o
+OBJS_SPEED = root/obj/node.o root/obj/graph.o examples/obj/speed.o
+OBJS_SIMPLE = root/obj/node.o root/obj/graph.o examples/obj/simple.o
+
+all : gradient_descent ann gradient speed simple
+
+gradient_descent : $(BIN) root/include/vectmath.h root/include/mor.h root/include/dor.h root/include/por.h
+	$(MAKE) -C examples obj obj/gradient_descent.o
+	$(MAKE) -C root obj obj/node.o obj/graph.o
+	$(CXX) -o $(BIN)/gradient_descent $(OBJS_GRADIENT_DESCENT) $(LIBS)
+
+ann : $(BIN) root/include/vectmath.h root/include/mor.h root/include/dor.h root/include/por.h
+	$(MAKE) -C examples obj obj/ann.o
+	$(MAKE) -C root obj obj/node.o obj/graph.o
+	$(CXX) -o $(BIN)/ann $(OBJS_ANN) $(LIBS)
+
+gradient : $(BIN) root/include/vectmath.h root/include/mor.h root/include/dor.h root/include/por.h
+	$(MAKE) -C examples obj obj/gradient.o
+	$(MAKE) -C root obj obj/node.o obj/graph.o
+	$(CXX) -o $(BIN)/gradient $(OBJS_GRADIENT) $(LIBS)
+
+speed : $(BIN) root/include/vectmath.h root/include/mor.h root/include/dor.h root/include/por.h
+	$(MAKE) -C examples obj obj/speed.o
+	$(MAKE) -C root obj obj/node.o obj/graph.o
+	$(CXX) -o $(BIN)/speed $(OBJS_SPEED) $(LIBS)
+
+simple : $(BIN) root/include/mor.h root/include/dor.h root/include/por.h
+	$(MAKE) -C examples obj obj/simple.o
+	$(MAKE) -C root obj obj/node.o obj/graph.o
+	$(CXX) -o $(BIN)/simple $(OBJS_SIMPLE) $(LIBS)
+
+$(BIN) :
+	if [ ! -d $(BIN) ]; then mkdir $(BIN); fi
+
+$(LIB) :
+	if [ ! -d $(LIB) ]; then mkdir $(LIB); fi
+
+clean :
+	$(MAKE) -C root clean
+	$(MAKE) -C examples clean
+	if [ -d $(BIN) ]; then rm $(BIN) -r; fi
+
+install : $(LIB)
+	ar rcs $(LIB)/lib$(LIB_NAME).a root/obj/graph.o root/obj/node.o
+	if [ ! -d /usr/local/include/$(LIB_NAME) ]; then sudo mkdir /usr/local/include/$(LIB_NAME); fi
+	sudo cp $(LIB)/lib$(LIB_NAME).a /usr/local/lib
+	sudo cp root/include/*.h /usr/local/include/$(LIB_NAME)
+
+.PHONY : all
+.PHONY : gradient_descent
+.PHONY : ann
+.PHONY : gradient
+.PHONY : speed
+.PHONY : simple
+.PHONY : clean
+.PHONY : install

install_diff/examples/Makefile

+CXX = g++
+ODIR = obj
+CXXFLAGS = -std=c++11 -O3
+OBJS = $(ODIR)/gradient_descent.o $(ODIR)/ann.o $(ODIR)/gradient.o $(ODIR)/speed.o $(ODIR)/simple.o 
+
+all : $(ODIR) $(OBJS) 
+
+$(ODIR)/gradient_descent.o : src/gradient_descent.cpp ../root/include/vectmath.h ../root/include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR)/ann.o : src/ann.cpp ../root/include/vectmath.h ../root/include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR)/gradient.o : src/gradient.cpp ../root/include/vectmath.h ../root/include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR)/speed.o : src/speed.cpp ../root/include/vectmath.h ../root/include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR)/simple.o : src/simple.cpp ../root/include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR) : 
+	if [ ! -d $(ODIR) ]; then mkdir $(ODIR); fi
+
+clean : 
+	if [ -d $(ODIR) ]; then rm $(ODIR) -r; fi
+
+.PHONY : all
+.PHONY : clean
+

install_diff/examples/src/ann.cpp

+#include <iostream>
+#include <ctime>
+
+#include "../../root/include/vectmath.h"
+#include "../../root/include/node.h"
+
+typedef std::vector<Node> Vector;
+typedef std::vector<Vector> Matrix;
+
+Node random_number(){
+    return rand()/(double)RAND_MAX;
+}
+
+Node tan_h(Node& x){
+    return (1-exp(-2*x))/(1+exp(-2*x));
+}
+
+Node mean_square_error(Vector& y_true, Vector& y_pred){
+    Node loss;
+    for(size_t i=0 ; i<y_true.size() ; i++){
+        loss += pow(y_true[i]-y_pred[i], 2);
+    }
+    return loss;
+}
+
+struct Layer {
+    Matrix weights;
+    Matrix bias;
+    Node (*activation)(Node&);
+    int input_shape;
+    int output_shape;
+
+    Layer(int input, int output, Node (*activation)(Node&)){
+        this->activation = activation;
+        this->input_shape = input;
+        this->output_shape = output;
+        weights.resize(input, Vector(output));
+        bias.resize(1, Vector(output));
+        random_number >> weights;
+        random_number >> bias;
+    }
+
+    Matrix forward(Matrix& previous){
+        Matrix output = dot(previous, weights) + bias;
+        output = activation >> output;
+        return output;
+    }
+
+    void backward(Node& loss, const float& learning_rate){
+        weights -= learning_rate*loss.gradient(weights);
+        bias -= learning_rate*loss.gradient(bias);
+    }
+};
+
+struct Network {
+    std::vector<Layer> layers;
+    int input_shape;
+    Graph* graph;
+
+    Network(){
+        graph = Graph::getInstance();
+    }
+
+    void input_layer(int input_shape){
+        this->input_shape = input_shape;
+    }
+
+    void add(int output_shape, Node (*activation)(Node&)){
+        int input = layers.empty()?input_shape:layers.back().output_shape;
+        layers.push_back(Layer(input, output_shape, activation));
+    }
+
+    Matrix run(Matrix& input){
+        Matrix output(input.size());
+        for(size_t j=0 ; j<input.size() ; j++){
+            Matrix out = {input[j]};
+            for(auto& lay : layers){
+                out = lay.forward(out);
+            }
+            output[j] = out[0];
+        }
+        return output;
+    }
+
+    void fit(Matrix& input, Matrix& output, Node (*loss_function)(Vector&, Vector&), int epochs, float learning_rate){
+        int p=0;
+        for(size_t i=0 ; i<epochs ; i++){
+            std::cout << "\r" << i+1 << "/" << epochs;
+            for(size_t j=0 ; j<input.size() ; j++){
+                // compute input
+                Matrix out = {input[j]};
+                for(auto& lay : layers){
+                    out = lay.forward(out);
+                }
+
+                // compute loss
+                Node loss = loss_function(output[j], out[0]);
+
+                // update parameters
+                for(auto& lay : layers){
+                    lay.backward(loss, learning_rate);
+                }
+
+                graph->new_recording();
+            }
+        }
+        std::cout << std::endl;
+    }
+};
+
+int main(int argc, char const *argv[]) {
+    srand(time(NULL));
+
+    Matrix input = {{0,0},{0,1},{1,0},{1,1}};
+    Matrix output = {{0},{1},{1},{0}};
+
+    Network network;
+    network.input_layer(2);
+    network.add(3, tan_h);
+    network.add(1, tan_h);
+    network.fit(input, output, mean_square_error, 500, 0.1);
+
+    Matrix pred = network.run(input);
+    std::cout << pred << std::endl;
+    return 0;
+}

install_diff/examples/src/gradient.cpp

+#include <iostream>
+
+#include "../../root/include/vectmath.h"
+#include "../../root/include/node.h"
+
+Node function(std::vector<Node> x){
+    return pow(x[0]-x[1], 2) + x[0]*x[1]*x[2]; // (x-y)^2 + x*y*z
+}
+
+int main(int argc, char const *argv[]) {
+    std::vector<Node> x = {5,6,7};
+    Node f = function(x);
+    std::cout << "grad(f) = " << f.gradient(x) << std::endl;
+    return 0;
+}

install_diff/examples/src/gradient_descent.cpp

+#include <iostream>
+
+#include "../../root/include/vectmath.h"
+#include "../../root/include/node.h"
+
+Node function(std::vector<Node>& x){
+    return pow(x[0], 2) + pow(x[1], 2); // x^2 + y^2
+}
+
+int main(int argc, char const *argv[]) {
+    Graph* graph = Graph::getInstance();
+
+    std::vector<Node> x = {50, 50};
+    Node f;
+
+    int epochs = 30;
+    float learning_rate = 0.1;
+    for(size_t i=0 ; i<epochs ; i++){
+        f = function(x);
+        x -= learning_rate*f.gradient(x);
+        graph->new_recording();
+    }
+
+    std::cout << "f = " << f << std::endl;
+    std::cout << "x = " << x << std::endl;
+    return 0;
+}

install_diff/examples/src/simple.cpp

+#include <iostream>
+#include "../../root/include/node.h"
+
+int main(int argc, char const *argv[]) {
+    Node x=2, y=3;
+    Node f = x*y + sin(x);
+
+    std::cout << "f(x,y) = x*y + sin(x)" << std::endl;
+    std::cout << "f(" << x << "," << y << ") = " << f << std::endl;
+    std::cout << "∂f/∂x = " << f.gradient(x) << std::endl;
+    std::cout << "∂f/∂y = " << f.gradient(y) << std::endl;
+    return 0;
+}

install_diff/examples/src/speed.cpp

+#include <iostream>
+#include <iomanip>
+#include <cassert>
+#include <ctime>
+#include <chrono>
+
+#include "../../root/include/vectmath.h"
+#include "../../root/include/node.h"
+
+template <class T>
+std::vector<std::vector<T> > get_random_matrix(const int& height, const int& width, T t){
+    std::vector<std::vector<T> > mat(height, std::vector<T>(width));
+    for(auto& v : mat){
+        for(auto& e : v){
+            e = rand()/(double)RAND_MAX;
+        }
+    }
+    return mat;
+}
+
+int main(int argc, char const *argv[]) {
+    srand(time(0));
+
+    int size = 30;
+    std::vector<std::vector<double> > a = get_random_matrix(size, size, double());
+    std::vector<std::vector<double> > b = get_random_matrix(size, size, double());
+    std::vector<std::vector<Node> > c = get_random_matrix(size, size, Node());
+    std::vector<std::vector<Node> > d = get_random_matrix(size, size, Node());
+
+    std::cout << std::fixed;
+    std::cout << std::setprecision(10);
+
+    std::cout << "Running with double...\t";
+    std::cout.flush();
+    auto start = std::chrono::high_resolution_clock::now();
+    std::vector<std::vector<double> > ab = dot(a, b);
+    auto finish = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> elapsed = finish - start;
+    std::cout << "Elapsed time: " << elapsed.count() << " s" << std::endl;
+
+    std::cout << "Running with Node...\t";
+    std::cout.flush();
+    start = std::chrono::high_resolution_clock::now();
+    std::vector<std::vector<Node> > cd = dot(c, d);
+    finish = std::chrono::high_resolution_clock::now();
+    elapsed = finish - start;
+    std::cout << "Elapsed time: " << elapsed.count() << " s" << std::endl;
+
+    std::cout << "Yet to be improved..." << std::endl;
+    return 0;
+}

install_diff/root/Makefile

+CXX = g++
+ODIR = obj
+CXXFLAGS = -std=c++11 -O3
+OBJS = $(ODIR)/graph.o $(ODIR)/node.o 
+
+all : $(ODIR) $(OBJS) 
+
+$(ODIR)/graph.o : src/graph.cpp include/graph.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR)/node.o : src/node.cpp include/node.h 
+	$(CXX) -c $< -o $@ $(CXXFLAGS)
+
+$(ODIR) : 
+	if [ ! -d $(ODIR) ]; then mkdir $(ODIR); fi
+
+clean : 
+	if [ -d $(ODIR) ]; then rm $(ODIR) -r; fi
+
+.PHONY : all
+.PHONY : clean
+

install_diff/root/include/dor.h

+#ifndef DYADIC_OPERATION_RESULT
+#define DYADIC_OPERATION_RESULT
+
+struct DyadicOperationResult {
+    double value;
+    double left_grad;
+    double right_grad;
+
+    DyadicOperationResult(double value, double left_grad, double right_grad){
+        this->value = value;
+        this->left_grad = left_grad;
+        this->right_grad = right_grad;
+    }
+};
+
+#endif /* end of include guard: DYADIC_OPERATION_RESULT */

install_diff/root/include/graph.h

+#ifndef GRAPH_H
+#define GRAPH_H
+
+#include <map>
+#include <vector>
+#include <utility>
+
+class Graph {
+    private:
+        std::map<long int, std::vector<std::pair<double, long int> > > nodes;
+        static Graph* instance;
+        Graph();
+
+    public:
+        static long int uid_counter;
+        static long int uid();
+        static Graph* getInstance();
+
+        void connect(const long int& uid, const std::pair<double, long int>& edge);
+        std::vector<std::pair<double, long int> > get(const long int& uid) const;
+        bool has(const long int& uid) const;
+
+        void new_recording();
+};
+
+#endif /* end of include guard: GRAPH_H */

install_diff/root/include/mor.h

+#ifndef MONADIC_OPERATION_RESULT
+#define MONADIC_OPERATION_RESULT
+
+struct MonadicOperationResult {
+    double value;
+    double grad;
+
+    MonadicOperationResult(double value, double grad){
+        this->value = value;
+        this->grad = grad;
+    }
+};
+
+#endif /* end of include guard: MONADIC_OPERATION_RESULT */

install_diff/root/include/node.h

+#ifndef NODE_H
+#define NODE_H
+
+#include <cmath>
+#include <iostream>
+
+#include "graph.h"
+#include "mor.h"
+#include "dor.h"
+#include "por.h"
+
+class Node {
+    private:
+        double value;
+        long int uid;
+
+        double gradient_recursive(Graph* graph, const long int& current_uid, const long int& stop_uid) const;
+
+    public:
+        Node(const double& value=0);
+        Node(const Node& node);
+
+        static Node monadic_operation(const Node& n, MonadicOperationResult (*)(const double&));
+        static Node dyadic_operation(const Node& l, const Node& r, DyadicOperationResult (*)(const double&, const double&));
+        static Node polyadic_operation(const std::vector<Node>& nodes, PolyadicOperationResult (*)(const std::vector<double>&));
+
+        double gradient(const Node& node) const;
+        std::vector<double> gradient(const std::vector<Node>& nodes) const;
+        std::vector<std::vector<double> > gradient(const std::vector<std::vector<Node> >& nodes) const;
+
+        friend Node operator+(const Node& l, const Node& r);
+        friend Node operator-(const Node& l, const Node& r);
+        friend Node operator*(const Node& l, const Node& r);
+        friend Node operator/(const Node& l, const Node& r);
+
+        Node& operator+=(const Node& r);
+        Node& operator-=(const Node& r);
+        Node& operator*=(const Node& r);
+        Node& operator/=(const Node& r);
+
+        friend bool operator==(const Node& l, const Node& r);
+        friend bool operator<(const Node& l, const Node& r);
+        friend bool operator>(const Node& l, const Node& r);
+        friend bool operator<=(const Node& l, const Node& r);
+        friend bool operator>=(const Node& l, const Node& r);
+
+        friend Node sin(const Node& x);
+        friend Node cos(const Node& x);
+        friend Node tan(const Node& x);
+        friend Node sinh(const Node& x);
+        friend Node cosh(const Node& x);
+        friend Node tanh(const Node& x);
+        friend Node asin(const Node& x);
+        friend Node acos(const Node& x);
+        friend Node atan(const Node& x);
+
+        friend Node log(const Node& x, const Node& base);
+        friend Node log10(const Node& x);
+        friend Node ln(const Node& x);
+
+        friend Node pow(const Node& x, const Node& p);
+        friend Node exp(const Node& x);
+        friend Node sqrt(const Node& x);
+
+        friend Node abs(const Node& x);
+        friend Node min(const Node& l, const Node& r);
+        friend Node max(const Node& l, const Node& r);
+
+        friend std::ostream& operator<<(std::ostream& os, const Node& node);
+};
+
+#endif /* end of include guard: NODE_H */

install_diff/root/include/por.h

+#ifndef POLYADIC_OPERATION_RESULT
+#define POLYADIC_OPERATION_RESULT
+
+#include <vector>
+
+struct PolyadicOperationResult {
+    double value;
+    std::vector<double> gradients;
+
+    PolyadicOperationResult(double value, const std::vector<double>& gradients){
+        this->value = value;
+        this->gradients = gradients;
+    }
+};
+
+#endif /* end of include guard: POLYADIC_OPERATION_RESULT */

install_diff/root/include/vectmath.h

+#ifndef VECTMATH
+#define VECTMATH
+
+#include <vector>
+#include <cassert>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+
+// dot product
+template <typename T>
+std::vector<std::vector<T> > dot(const std::vector<std::vector<T> >& a, const std::vector<std::vector<T> >& b){
+    assert(a[0].size()==b.size());
+
+    T w=0;
+    std::vector<std::vector<T> > result(a.size(), std::vector<T>(b[0].size()));
+    for (int i=0 ; i<a.size() ; i++){
+        for (int j=0 ; j<b[0].size() ; j++){
+            for (int h=0 ; h<b.size() ; h++){
+                w += a[i][h]*b[h][j];
+            }
+            result[i][j] = w;
+            w=0;
+        }
+    }
+
+    return result;
+}
+
+// operators
+template <typename U, typename V>
+std::vector<U>& operator-=(std::vector<U>& u, const std::vector<V>& v){
+    assert(u.size()==v.size());
+    for(size_t i=0 ; i<u.size() ; i++){
+        u[i] -= v[i];
+    }
+    return u;
+}
+
+template <typename U>
+std::vector<U> operator+(const std::vector<U>& u, const std::vector<U>& v){
+    assert(u.size()==v.size());
+    std::vector<U> w(u.size());
+    for(size_t i=0 ; i<w.size() ; i++){
+        w[i] = u[i]+v[i];
+    }
+    return w;
+}
+
+template <typename U, typename S>
+std::vector<U> operator*(const S& s, const std::vector<U>& u){
+    std::vector<U> result(u.size());
+    for(size_t i=0 ; i<u.size() ; i++){
+        result[i] = s*u[i];
+    }
+    return result;
+}
+
+template <typename U>
+std::vector<U>& operator>>(U (*fun)(U&), std::vector<U>& u){
+    std::transform(u.begin(), u.end(), u.begin(), fun);
+    return u;
+}
+
+template <typename U, typename S>
+std::vector<U>& operator>>(S (*fun)(S&), std::vector<U>& u){
+    for(auto& v : u){
+        fun >> v;
+    }
+    return u;
+}
+
+template <typename U>
+std::vector<U>& operator>>(U (*fun)(), std::vector<U>& u){
+    for(auto& e : u){
+        e = fun();
+    }
+    return u;
+}
+
+template <typename U, typename S>
+std::vector<U>& operator>>(S (*fun)(), std::vector<U>& u){
+    for(auto& v : u){
+        fun >> v;
+    }
+    return u;
+}
+
+template <typename U>
+std::ostream& operator<<(std::ostream& os, const std::vector<U>& u){
+    os << "[";
+    for(size_t i=0 ; i<u.size() ; i++){
+        os << u[i] << (i<u.size()-1?", ":"");
+    }
+    os << "]";
+    return os;
+}
+
+#endif

install_diff/root/src/graph.cpp

+#include "../include/graph.h"
+
+Graph* Graph::instance = 0;
+long int Graph::uid_counter = 1;
+
+long int Graph::uid(){
+    return uid_counter++;
+}
+
+Graph* Graph::getInstance(){
+    if(instance==NULL){
+        instance = new Graph();
+    }
+    return instance;
+}
+
+Graph::Graph(){}
+
+void Graph::connect(const long int& uid, const std::pair<double, long int>& edge){
+    nodes[uid].push_back(edge);
+}
+
+std::vector<std::pair<double, long int> > Graph::get(const long int& uid) const{
+    return nodes.at(uid);
+}
+
+bool Graph::has(const long int& uid) const{
+    return nodes.find(uid)!=nodes.end();
+}
+
+void Graph::new_recording(){
+    nodes.clear();
+}

install_diff/root/src/node.cpp

+#include "../include/node.h"
+
+Node::Node(const double& value) {
+    this->value = value;
+    this->uid = Graph::uid();
+}
+
+Node::Node(const Node& node){
+    this->value = node.value;
+    this->uid = node.uid;
+}
+
+double Node::gradient_recursive(Graph* graph, const long int& current_uid, const long int& stop_uid) const{
+    if(current_uid==stop_uid){
+        return 1.0;
+    }
+
+    double sum=0.0;
+    if(graph->has(current_uid)){
+        for(auto& pair : graph->get(current_uid)){
+            sum += pair.first*gradient_recursive(graph, pair.second, stop_uid);
+        }
+    }
+    return sum;
+}
+
+double Node::gradient(const Node& node) const{
+    Graph* graph = Graph::getInstance();
+    return gradient_recursive(graph, this->uid, node.uid);
+}
+
+std::vector<double> Node::gradient(const std::vector<Node>& nodes) const{
+    Graph* graph = Graph::getInstance();
+    std::vector<double> grad(nodes.size());
+    for(size_t i=0 ; i<nodes.size() ; i++){
+        grad[i] = gradient_recursive(graph, this->uid, nodes[i].uid);
+    }
+    return grad;
+}
+
+std::vector<std::vector<double> > Node::gradient(const std::vector<std::vector<Node> >& nodes) const{
+    Graph* graph = Graph::getInstance();
+    std::vector<std::vector<double> > grad(nodes.size());
+    for(size_t i=0 ; i<nodes.size() ; i++){
+        grad[i].resize(nodes[i].size());
+        for(size_t j=0 ; j<nodes[i].size() ; j++){
+            grad[i][j] = gradient_recursive(graph, this->uid, nodes[i][j].uid);
+        }
+    }
+    return grad;
+}
+
+Node Node::monadic_operation(const Node& n, MonadicOperationResult (*fun)(const double&)){
+    MonadicOperationResult res = fun(n.value);
+    Node result(res.value);
+    Graph* graph = Graph::getInstance();
+    graph->connect(result.uid, std::make_pair(res.grad, n.uid));
+    return result;
+}
+
+Node Node::dyadic_operation(const Node& left, const Node& right, DyadicOperationResult (*fun)(const double&, const double&)){
+    DyadicOperationResult res = fun(left.value, right.value);
+    Node result(res.value);
+    Graph* graph = Graph::getInstance();
+    graph->connect(result.uid, std::make_pair(res.left_grad, left.uid));
+    graph->connect(result.uid, std::make_pair(res.right_grad, right.uid));
+    return result;
+}
+
+Node Node::polyadic_operation(const std::vector<Node>& nodes, PolyadicOperationResult (*fun)(const std::vector<double>&)){
+    std::vector<double> values(nodes.size());
+    for(size_t i=0 ; i<nodes.size() ; i++){
+        values[i] = nodes[i].value;
+    }
+    PolyadicOperationResult res = fun(values);
+    Node result(res.value);
+    Graph* graph = Graph::getInstance();
+    for(size_t i=0 ; i<nodes.size() ; i++){
+        graph->connect(result.uid, std::make_pair(res.gradients[i], nodes[i].uid));
+    }
+    return result;
+}
+
+Node operator+(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& l, const double& r){
+        return DyadicOperationResult(l+r, 1.0, 1.0);
+    });
+}
+
+Node operator-(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& l, const double& r){
+        return DyadicOperationResult(l-r, 1.0, -1.0);
+    });
+}
+
+Node operator*(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& l, const double& r){
+        return DyadicOperationResult(l*r, r, l);
+    });
+}
+
+Node operator/(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& l, const double& r){
+        return DyadicOperationResult(l/r, 1.0/r, -1.0*l/(r*r));
+    });
+}
+
+Node& Node::operator+=(const Node& r){
+    *this = *this + r;
+    return *this;
+}
+
+Node& Node::operator-=(const Node& r){
+    *this = *this - r;
+    return *this;
+}
+
+Node& Node::operator*=(const Node& r){
+    *this = *this * r;
+    return *this;
+}
+
+Node& Node::operator/=(const Node& r){
+    *this = *this / r;
+    return *this;
+}
+
+bool operator==(const Node& left, const Node& right){
+    return left.value==right.value;
+}
+
+bool operator<(const Node& left, const Node& right){
+    return left.value<right.value;
+}
+
+bool operator>(const Node& left, const Node& right){
+    return left.value>right.value;
+}
+
+bool operator<=(const Node& left, const Node& right){
+    return left.value<=right.value;
+}
+
+bool operator>=(const Node& left, const Node& right){
+    return left.value>=right.value;
+}
+
+Node sin(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::sin(n), ::cos(n));
+    });
+}
+
+Node cos(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::cos(n), -1.0*::sin(n));
+    });
+}
+
+Node tan(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::tan(n), 1.0/::pow(::cos(n), 2));
+    });
+}
+
+Node sinh(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::sinh(n), ::cosh(n));
+    });
+}
+
+Node cosh(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::cosh(n), ::sinh(n));
+    });
+}
+
+Node asin(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::asin(n), 1.0/(::sqrt(1-n*n)));
+    });
+}
+
+Node acos(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::acos(n), -1.0/(::sqrt(1-n*n)));
+    });
+}
+
+Node atan(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::atan(n), 1.0/(1+n*n));
+    });
+}
+
+Node tanh(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::tanh(n), 1.0-::pow(::tanh(n), 2));
+    });
+}
+
+Node log(const Node& x, const Node& base){
+    return Node::dyadic_operation(x, base, [](const double& a, const double& b){
+        return DyadicOperationResult(::log(a)/::log(b), 1.0/(a*::log(b)), -1.0*::log(a)/(b*::log(b)));
+    });
+}
+
+Node log10(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::log(n)/::log(10), 1.0/(n*::log(10)));
+    });
+}
+
+Node ln(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::log(n), 1.0/::log(n));
+    });
+}
+
+Node pow(const Node& x, const Node& base){
+    return Node::dyadic_operation(x, base, [](const double& a, const double& b){
+        if(a<=0){
+            return DyadicOperationResult(::pow(a,b), b*::pow(a,b-1), 0);
+        }
+        return DyadicOperationResult(::pow(a,b), b*::pow(a,b-1), ::log(a)*::pow(a,b));
+    });
+}
+
+Node exp(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::exp(n), ::exp(n));
+    });
+}
+
+Node sqrt(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        return MonadicOperationResult(::sqrt(n), 1.0/(2*::sqrt(n)));
+    });
+}
+
+Node abs(const Node& x){
+    return Node::monadic_operation(x, [](const double& n){
+        int sign = n==0 ? 0 : n/::abs(n);
+        return MonadicOperationResult(::abs(n), sign);
+    });
+}
+
+Node min(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& a, const double& b){
+        if(a<b){
+            return DyadicOperationResult(a, 1, 0);
+        }
+        if(a>b){
+            return DyadicOperationResult(b, 0, 1);
+        }
+        return DyadicOperationResult(a, 0, 0);
+    });
+}
+
+Node max(const Node& left, const Node& right){
+    return Node::dyadic_operation(left, right, [](const double& a, const double& b){
+        if(a>b){
+            return DyadicOperationResult(a, 1, 0);
+        }
+        if(a<b){
+            return DyadicOperationResult(b, 0, 1);
+        }
+        return DyadicOperationResult(a, 0, 0);
+    });
+}
+
+std::ostream& operator<<(std::ostream& os, const Node& node){
+    os << node.value;
+    return os;
+}

root_diff @ 77e0b021

-Subproject commit 77e0b021330e45393302d3b0301945230e708e50