AronLib.h

Go to the documentation of this file.

#ifndef __AronLib__
#define __AronLib__

#ifdef __cplus
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <utility> // std::pair, std::make_pair
#endif

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <utility> // std::pair, std::make_pair

#include <execinfo.h>
#include <stdio.h>
#include <stdlib.h>
#include "Vector3.h"
#include "DDLinkedList.h"
#include "Const.h"
#include <experimental/optional>
// #include "/Library/Developer/CommandLineTools/usr/include/c++/v1/experimental/optional"
#include <boost/filesystem/operations.hpp>

#include <vector>
#include </Users/cat/myfile/bitbucket/cpplib/glm/glm/glm.hpp>


#include <dirent.h>
#include <unistd.h>
#include <limits.h>
#include <regex>

#include <algorithm>
#include <cctype>
#include <locale>

#include <limits>    // for std::numeric_limits
#include <stdexcept> // for std::overflow_error

#include <array>  // std::array<int, 2> arr = {1, 2};

#ifdef __has_include                           // Check if __has_include is present
#  if __has_include(<optional>)                // Check for a standard library
#    include<optional>
#  elif __has_include(<experimental/optional>) // Check for an experimental version
#    include <experimental/optional>
#  elif __has_include(<boost/optional.hpp>)    // Try with an external library
#    include <boost/optional.hpp>
#  else                                        // Not found at all
#     error "Missing <optional>"
#  endif
#endif



//const string S_RED = "\033[1;31m";
//const string S_END = "\033[0m";

//cout << "\033[1;31mbold red text\033[0m\n";

namespace fs = boost::filesystem;

using namespace std;


template<typename T>
void pl(T t);

void nl();

template<typename T>
void swap(vector<T>& vec, int i, int j);

int random(int n);

void nl(){
    printf("\n");
}

int random(int n) {
    return arc4random_uniform(n);
}

// print all prime from 2 to n
vector<int> allPrime(int n) {
    vector<int> vec;
    if(n >= 2) {
        vec.push_back(2);
        for(int p=3; p<n; p++) {
            bool isPrime = true;
            for(int i=0; i<vec.size() && isPrime; i++) {
                if(p % vec[i] == 0)
                    isPrime = false;
            }
            if(isPrime)
                vec.push_back(p);
        }
    }
    return vec;
}



const char* stringToCString(string s) {
    return s.c_str();
}


const char* s2c(string s) {
    return s.c_str();
}

const char* toCharPtr(string s) {
    return s2c(s);
}

template<typename T>
string toStr(const T& s) {
    std::ostringstream out;
    out<<s;
    return out.str();
}

template<typename T, typename U>
string toStr(std::pair<T, U> p){
    string ret = "";
    ret += "(";
    ret += toStr(p.first) + " , ";
    ret += toStr(p.second);
    ret += ")";
    return ret;
}

template<typename T, typename U>
string str(std::pair<T, U> p){
    return toStr(p);
}

template<typename T>
string str(const T& s) {
    return toStr(s);
}

template<typename T>
string toStr2(T& s) {
    std::ostringstream out;
    out<<s;
    return out.str();
}

template<typename T>
string str2(T& s) {
    return toStr2(s);
}



string trimLeft(string str) {
    string s = str;
    s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](auto& x) {
        return !std::isspace(x);
    }));
    return s;
}

string trimRight(string str) {
    string s = str;
    s.erase(std::find_if(s.rbegin(), s.rend(), [](auto& x) {
        return !std::isspace(x);
    }).base(), s.end());
    return s;
}

string trim(string str) {
    string s = str;
    return trimRight(trimLeft(s));
}

string cStringToString(const char* pt) {
    std::string s(pt);
    return s;
}

string c2s(const char* pt) {
    string s(pt);
    return s;
}

bool containStrRegex(string s, regex rx) {
    smatch m;
    regex_search(s, m, rx);
    return m.size() > 0;
}

bool containStr(string str, string pat) {
    smatch m;
    regex rx(pat);
    regex_search(str, m, rx);
    return m.size() > 0;
}

vector<string> splitStr(string s, string delim) {
    vector<string> vec;
    auto start = 0U;
    auto end = s.find(delim);
    while (end != std::string::npos) {
        vec.push_back(s.substr(start, end - start));
        start = end + delim.length();
        end = s.find(delim, start);
    }
    vec.push_back(s.substr(start, s.size() - start));
    return vec;
}


template<typename T>
std::pair<vector<T>, vector<T>> splitAt(int n, vector<T> vec) {
    vector<T> fstVec;
    vector<T> sndVec;
    for(int i=0; i<vec.size(); i++) {
        if(i < n)
            fstVec.push_back(vec[i]);
        else
            sndVec.push_back(vec[i]);
    }
    return std::make_pair(fstVec, sndVec);
}


template<typename Fun, typename T>
std::pair<vector<T>, vector<T>> breakVec(Fun f, vector<T> vec) {
    bool first = true;
    vector<T> fstVec;
    for(int i=0; i<vec.size(); i++) {
        if(f(vec[i])) {
            return splitAt(i, vec);
        }
    }
    return std::make_pair(fstVec, vec);
}

template<typename T>
vector<T> replaceVec(T oldVal, T newVal, vector<T> vec) {
    vector<T> retVec(vec);
    std::replace(retVec.begin(), retVec.end(), oldVal, newVal);
    return retVec;
}

template<typename T>
vector<T> insertAt(vector<T> vec, int pos, vector<T> insertVec) {
    assert(0 <= pos && pos <= vec.size());
    vector<T> v(vec);
    typename vector<T>::iterator it = v.begin();
    advance(it, pos);
    std::copy (insertVec.begin(),insertVec.end(),std::inserter(v,it));
    return v;
}



vector<std::string> splitStrRegex(const string& s, string rgxStr = "\\s+") {
    vector<std::string> vec;
    regex rgx (rgxStr);

    sregex_token_iterator iter(s.begin(), s.end(), rgx, -1);
    sregex_token_iterator end;

    while (iter != end) {
        vec.push_back(*iter);
        ++iter;
    }
    return vec;
}



bool isEmpty(string s) {
    return s.empty();
}

bool isEmpty(const char* cpt) {
    if(cpt != NULL) {
        return strlen(cpt) == 0;
    } else {
        perror("cpt is NULL");
        exit(EXIT_FAILURE);
    }
    return false;
}

double squareRoot(double a) {
    double x0 = a + 1;
    double x = x0;
    while(fabs(x*x - a) > 0.00000001) {
        x = x + (a - x*x)/(2*x);
    }
    return x;
}

double intToDouble(int n) {
    return static_cast<double>(n);
}

double nthRoot(double c, int n, double epsilon = 0.000001) {

    // f(x) = x^{1/2}
    auto ff = [](auto xx, auto nn, auto cc) {
        double nb = intToDouble(nn);
        return pow(xx, nb) - cc;
    };

    // The first derivative of f(x)
    auto df = [](auto xx, auto nn, auto cc) {
        double bb = intToDouble(nn - 1);
        return nn*pow(xx, bb);
    };

    double b = intToDouble(n);
    double x = c + 1;
    while(fabs(pow(x, b) - c) > epsilon) {
        x = x - ff(x, n, c)/df(x, n, c);
    }
    return x;
}


// KEY: char array to string, char arr to str
string charArrToString(const char* pt) {
    return cStringToString(pt);
}


// TODO: find a better way to do it
// replicate n copy of s
string replicate(int n, string s) {
    string str = "";
    for(int i=0; i<n; i++) {
        str += s;
    }
    return str;
}

// TODO: find a better way to do it
// repeating n copy of s
string repeat(int n, string s) {
    return replicate(n, s);
}

template<typename T>
vector<T> repeatVec(int n, T t) {
    vector<T> vec;
    for(int i=0; i<n; i++) {
        vec.push_back(t);
    }
    return vec;
}


int compareCString(char* s1, char* s2) {
    return strcmp(s1, s2);
}

int compareString(string s1, string s2) {
    return s1.compare(s2);
}

bool compareFloat(double a, double b) {
    double epsilon = 0.00001;
    return fabs(a-b) < epsilon;
}

void print(glm::vec3 v3) {
    printf("%2.3f %2.3f %2.3f \n",v3.x, v3.y, v3.z);
}
void print(glm::vec2 v2) {
    printf("%2.3f %2.3f \n",v2.x, v2.y);
}

void printLn(glm::vec3 v3) {
    printf("%2.3f %2.3f %2.3f ",v3.x, v3.y, v3.z);
}
void printLn(glm::vec2 v2) {
    printf("%2.3f %2.3f ",v2.x, v2.y);
}


string toString(glm::vec3 v3) {
    char strPt[50];
    sprintf(strPt, "%2.3f %2.3f %2.3f", v3.x, v3.y, v3.z);
    return charArrToString(strPt);
}


template<typename T>
vector<vector<T>> transpose(vector<vector<T>> mat) {
    vector<vector<T>> retMat;
    int min = INT_MAX;
    for(auto const& v : mat) {
        if(min > v.size())
            min = v.size();
    }

    for(int i=0; i<min; i++) {
        vector<T> v;
        for(int k = 0; k < mat.size(); k++) {
            v.push_back(mat[k][i]);
        }
        retMat.push_back(v);
    }
    return retMat;
}

// Compare array, compare int array
template<typename T>
bool compareArray(T arr1[], T arr2[], const int s1) {
    for(int i=0; i<s1; i++) {
        if(arr1[i] != arr2[i])
            return false;
    }
    return true;
}

template<typename T>
void print(std::pair<T, T> pair) {
    cout<<"("<<pair.first<<" "<<pair.second<<")"<<endl;
}

template<typename T>
void print(vector<T> vec) {
    for(auto const& v: vec) {
        std::cout<<"["<<v<<"]";
    }
    std::cout<<std::endl;
}

template<typename T>
void print(vector< vector<T> > vec) {
    for(auto const& v: vec) {
        print(v);
    }
    std::cout<<std::endl;
}

template<typename T>
void print(T t) {
    std::cout<<"["<<t<<"]";
}

template<typename T>
void printNb(T t) {
    std::cout<<t;
}

template<typename T>
void printLn(T t) {
    std::cout<<"["<<t<<"]"<<std::endl;
}

template<typename T>
void printLnNb(T t) {
    std::cout<<t<<std::endl;
}


template<typename T>
void pl(T t) {
    print(t);
    std::cout<<std::endl;
}

template<typename T>
void ppL(T t) {
    printLn(t);
}



template<typename T>
vector<T> operator+(vector<T> v1, vector<T> v2) {
    vector<T> vec(v1);
    vec.insert(vec.end(), v2.begin(), v2.end());
    return vec;
}

string operator+(const char* chars, string s) {
    return charArrToString(chars) + s;
}

string operator+(string s, const char* chars) {
    return s + charArrToString(chars);
}

template<typename T>
vector<T> con(T t, vector<T> v) {
    vector<T> vec;
    vec.push_back(t);
    for(auto const& e: v)
        vec.push_back(e);
    return vec;
}


template<typename T>
vector<T>* moveVec(T a, vector<T> vec) {
    vector<T>* pt = new vector<T>();
    for(auto const& val: vec) {
        pt -> push_back(val);
    }
    pt -> push_back(a);
    return std::move(pt);
}

int unsignToInt(size_t u){
    if(u > std::numeric_limits<int>::max()){
        throw std::overflow_error("ERROR: size_t > int");
    }
    return static_cast<int>(u);
}

template<typename T>
vector<T>* append(vector<T> vec, T a) {
    vector<T>* pt = new vector<T>();
    for(auto const& val: vec) {
        pt -> push_back(val);
    }
    pt -> push_back(a);
    return std::move(pt);
}

int length(string s){
    return unsignToInt(s.size());
}


int len(string s){
    return length(s);
}

template<typename T>
int length(vector<T> v){
    return unsignToInt(v.size());
}

template<typename T>
int len(vector<T> v){
    return length(v);
}

template<typename T>
bool binSearch(T key, T arr[], int lo, int hi){
    if(lo < hi){
        int mid = (lo + hi)/2;
        if(key < arr[mid]){
            return binSearch(key, arr, lo, mid - 1);    
        }else{
            return binSearch(key, arr, mid + 1, hi);
        }
    }else if(lo == hi){
        return key == arr[lo];
    }
    return false;
}

template<typename T>
vector<T> operator*(T a, vector<T> vec) {
    vector<T> newVec;
    for(auto const& value: vec) {
        newVec.push_back(a*value);
    }
    return newVec;
}

template<typename T>
vector<T> operator*(vector<T> vec, T a) {
    vector<T> newVec;
    for(auto const& value: vec) {
        newVec.push_back(a*value);
    }
    return newVec;
}

template<typename T>
vector<vector<T>> operator*(T a, vector<vector<T>> vec2) {
    vector<vector<T>> newVec2;
    vector<T> vec;
    for(auto const& vec : vec2) {
        newVec2.push_back(a*vec);
    }
    return newVec2;
}

template<typename T>
vector<vector<T>> operator*(vector<vector<T>> vec2, T a) {
    return a*vec2;
}




// Thanks for the default parameter
vector<int> range(int a, int b, int stride = 1) {
    vector<int> vec;
    for(int i=a; i<=b; i += stride)
        vec.push_back(i);
    return vec;
}

// generate vector double
template<typename T>
vector<T> geneVector(int a, int b) {
    int len = b - a;
    std::vector<T> vec;
    for(int i=0; i<=len; i++) {
        vec.push_back(a + i);
    }
    return vec;
}

template<typename T>
vector< vector<T> > geneMatrix(int ncol, int nrow, T init) {
    vector< vector<T> > vv;
    T val = init;
    for(int c=0; c<ncol; c++) {
        vector<T> tv;
        for(int r=0; r<nrow; r++) {
            tv.push_back(val);
            val++;
        }
        vv.push_back(tv);
    }
    return vv;
}


// swap array
template<typename T>
void swapT(T array[], int i, int j) {
    T tmp = array[i];
    array[i] = array[j];
    array[j] = tmp;
}

template<typename T>
T partitionT(T array[], int lo, int hi) {
    if(hi > lo) {
        T p = array[hi];
        T big = lo;
        for(int i=lo; i<=hi; i++) {
            if(array[i] <= p)     {
                swapT(array, i, big);
                if(i < hi)
                    big++;
            }
        }
        return big;
    }
    return -1;
}


// 2 3
// i = 0
//   j=1

template<typename T>
T partition(vector<T>& vec, int lo, int hi){
    if(lo < hi){
        T pivot = vec[hi];
        vector<T> tmp(hi - lo + 1);
        int i = lo;
        int j = hi;
        for(int k = 0; k < len(vec) && i != j; k++){
            if(vec[k] < pivot){
                tmp[i] = vec[k];
                i++;
            }else{
                tmp[j] = vec[k];
                j--;
            }
        }
        tmp[i] = pivot;
                
        for(int k=0; k < len(tmp); k++)
            vec[lo + k] = tmp[k];
        
        return i;
    }else{
        return lo;
    }
}


template<typename T>
T partititonInline(vector<T>& vec, int lo, int hi){
    int bigInx = lo;
    if(lo < hi){
        T pivot = vec[hi];
        for(int i = lo; i <= hi; i++){
            if(vec[i] <= pivot){
                swap(vec, i, bigInx);
                if(i != hi)
                    bigInx++;
            }else if(vec[i] > pivot){
                // vec[i] >= pivot
            }
        }
        
    }
    return bigInx;
}


// squick sort
template<typename T>
void quickSortT(T array[], int lo, int hi) {
    if(hi > lo) {
        T pivot = partitionT(array, lo, hi);
        quickSortT(array, lo, pivot-1);
        quickSortT(array, pivot+1, hi);
    }
}


class StopWatch {
    clock_t startTime = 0;
    clock_t endTime = 0;
    clock_t diff = 0;
public:
    StopWatch() {
        startTime = clock();
    }

    // Deprecated, use start() instead
    void begin() {
        startTime = clock();
        pl("Start StopWatch()");
    }
    void start() {
        begin();
    }
    void end() {
        endTime = clock();
        diff = endTime - startTime;
        std::cout<<"Diff=["<<(float)diff/CLOCKS_PER_SEC<<"] Seconds"<<std::endl;
        std::cout<<"Tick=["<<diff<<"] ticks"<<std::endl;
    }
    void print() {
        endTime = clock();
        diff = endTime - startTime;
        std::cout<<"Diff=["<<(float)diff/CLOCKS_PER_SEC<<"] Seconds"<<std::endl;
        std::cout<<"Tick=["<<diff<<"] ticks"<<std::endl;
    }
    void printTicket() {
        endTime = clock();
        clock_t diff = endTime - startTime;
        std::cout<<"Tick=["<<diff<<"] CPU tickets"<<std::endl;
    }
};


namespace SpacePrint {
    void pp(const char* format, ...);
    void fl();
    void fw(string s);

    template<typename T>
    void pp(T t);

    template<typename T>
    void pp(T t1, T t2);

    template<typename T>
    void pp(T t1, T t2, T t3);

    template<typename T, typename U>
    void pp(std::pair<T, U> p); 

    void pp(int n); 
    void pp(float n);
    void pp(long n); 
    void pp(unsigned long n); 


    template<typename T>
    void pp(T t) {
        print(t);
    }

    template<typename T>
    void pp(T t1, T t2) {
        print(t1);
        print(t2);
    }

    template<typename T>
    void pp(T t1, T t2, T t3) {
        print(t1);
        print(t2);
        print(t3);
    }

    template<typename T, typename U>
    void pp(std::pair<T, U> p){ 
        pp(toStr(p));
    }

    void pp(const char* format, ...) {
        printf(format);
    }

    void pp(int n) {
        printf("[%d]", n);
    }

    void pp(float n) {
        printf("[%f6]", n);
    }

    void pp(long n) {
        printf("[%ld]", n);
    }
    void pp(unsigned long n) {
        printf("[%lu]", n);
    }

    

    void fl() {
        pl("");
        printLnNb(repeat(80, "-"));
    }
    void fl(string s) {
        int n = 80;
        string s1 = repeat(n/2 - len(s)/2, "-");
        string s2 = repeat(n/2 - len(s)/2, "-");
        printLnNb(s1 + s + s2);
    }
    void fw(string s) {
        int len = s.length(); 
        int pad = 40 - (len/2);
        printLnNb(repeat(pad, "-") + s + repeat(pad, "-"));
    }

} // end SpacePrint

namespace SpaceTest {
using namespace SpacePrint;

    static void t(double a, double b, string s = "") {
        bool ret = compareFloat(a, b);
        if(ret)
            pp("true => " + s);
        else{
            pp(a, b);
            pp("false => " + s);
        }
        nl();
    }
    static void t(int a, int b, string s = "") {
       if(a == b) 
           pp("true => " + s);
       else{
           pp(a, b);
           pp("false => " + s);
       }
       nl();
    }
    static void t(long a, long b, string s = "") {
       if(a == b) 
           pp("true => " + s);
       else{
           pp(a, b);
           pp("false => " + s);
       }
       nl();
    }
    static void t(bool a, bool b, string s = "") {
       if(a == b) 
           pp("ture => " + s);
       else{
           pp(a, b);
           pp("false => " + s);
       }
       nl();
    }
    static void t(string a, string b, string s = "") {
       if(a == b) 
           pp("true => " + s);
       else{
           pp(a, b);
           pp("false => " + s);
       }
       nl();
    }
};

template<typename T>
void swap(T arr[], int i, int j) {
    T tmp = arr[i];
    arr[i] = arr[j];
    arr[j] = tmp;
}

template<typename T>
void swap(vector<T>& vec, int i, int j){
    T tmp = vec[i];
    vec[i] = vec[j];
    vec[j] = tmp;
}

// all the geometry functions here
namespace AronGeometry {
using namespace SpacePrint;

template<class T = double> class Vector;

template<class T = double>
class Point {
public:
    T x;
    T y;
    T z;
    Point() {
        x = 0.0;
        y = 0.0;
        z = 0.0;
    }
    Point(T x, T y, T z = 0.0 ) {
        this->x = x;
        this->y = y;
        this->z = z;
    }

    friend Point<T> operator+(Point<T> p1, Point<T> p2) {
        Point<T> p(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z);
        return p;
    }
    friend Point<T> operator-(Point<T> p1, Point<T> p2) {
        Point<T> p(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z);
        return p;
    }
    friend Point<T> operator*(double k, Point<T> p) {
        Point<T> p1(k*p.x, k*p.y, k*p.z);
        return p1;
    }
    friend Point<T> operator*(Point<T> p, double k) {
        return k*p;
    }


//        /*
//            vector(p0, p1) = p1 - p0
//        */
//        // TODO: add test
//        Vector<T> operator-(Point<T> p1){
//            Point<T> zero(0, 0, 0);
//            Point<T> p(p1.x - x, p1.y - y, p1.z - z);
//            Vector<T> v(zero, p);
//            return v;
//        }


    /*
        p1 = p0 + (p1 - p0) = p0 + v
    */
    // TODO: add test
    Point<T> operator+(Vector<T> v) {
        Point p(x + v.x, y + v.y, z + v.z);
        return p;
    }
    bool operator==(Point p0) {
        return x == p0.x && y == p0.y && z == p0.z;
    }
    bool operator!=(Point p0) {
        return x != p0.x || y != p0.y || z != p0.z;
    }
};


/*
    /Library/WebServer/Documents/zsurface/pdf/colinear_point.pdf

    p on segment p1 p2 => return 0
    p on one side of segment p1 p2 => return > 0
    p on other size of segment p1 p2 => return < 0
*/

template<class T>
class Pair {
public:
    Point<T> p0;
    Point<T> p1;
    Pair(Point<T> p0, Point<T> p1) {
        this -> p0 = p0;
        this -> p1 = p1;
    }
public:
    Pair operator+(Pair pair) {
        Pair pa(this -> p0 + pair.p0, this -> p1 + pair.p1);
        return pa;
    }
    Pair operator-(Pair pair) {
        Pair pa(this -> p0 - pair.p0, this -> p1 - pair.p1);
        return pa;
    }
};



template<class T>
class Vector {
public:
    Point<T> p0; // origin
    Point<T> p1;
public:
    Vector(Point<T> p0, Point<T> p1) {
        this -> p0 = p0;
        this -> p1 = p1;
    }

    bool isZero() {
        return p0 == p1;
    }
    Vector<T> dir() {
        // affine formula
        // vector: p0 -> p1 = p1 - p0
        return p1 - p0;
    }

    // v1 + v1' =
    // v1 = p1 - p0
    // v1' = p1' - p0'
    // v1 + v1' = (p1 - p0) + (p1' - p0')
    // Affine axiom:
    // there is unique map f such as
    // f:(S, V) -> S
    // a + v = bf
    Vector<T> operator+(Vector<T> v) {
        Vector<T> v1(p0 + v.p0, p1, v.p1);
        return v1;
    }
};

template<class T = double>
double isColinear3(Point<T> p, Point<T> p1, Point<T> p2) {
    // TODO: ignore z for now
    return (p2.y - p.y)*(p2.x - p1.x) - (p2.x - p.x)*(p2.y - p1.y);
}
template<class T = double>
bool isCrossSegment(Point<T> p0, Point<T> p1) {
    bool ret = false;
    return ret;
}

template<class T = double>
class Segment {
public:
    Point<T> p0;
    Point<T> p1;
public:
    /*
        Segment definition:
        1. p0 != p2
    */
    Segment(Point<T> p0, Point<T> p1) {
        if(p0 != p1) {
            this -> p0 = p0;
            this -> p1 = p1;
        } else {
            pp("some error");
        }
    }

    double squareDist() {
        double d = (p1.x - p0.x)*(p1.x - p0.x) +
                   (p1.y - p0.y)*(p1.y - p0.y);
        return d;
    }
    double distance() {
        double d = (p1.x - p0.x)*(p1.x - p0.x) +
                   (p1.y - p0.y)*(p1.y - p0.y);
        return sqrt(d);
    }
    bool isColinear(Point<T> p) {
        return isColinear3(p, p0, p1) == 0;
    }

    /*
        1. If two segment are the same, return false
        2. If one end point of segment on other segment, return true
        3.
    */
    bool isCrossSegment(Segment<T> s1, Segment<T> s2) {
        return true;
    }
    bool isEndPoint(Point<T> p) {
        return p == p0 || p == p1;
    }
    bool operator==(Segment<T> s) {
        return (p0 == s.p0 && p1 = s.p1) || (p1 = s.p0 && p0 == s.p1);
    }
};

}



void writeFile(string fname, string s) {
    ofstream ofs;
    ofs.open(fname, std::ofstream::out);
    ofs<<s;
    ofs.close();
}
void writeFileAppend(string fname, string s) {
    ofstream ofs;
    ofs.open(fname, std::ofstream::out | std::ofstream::app);
    ofs<<s;
    ofs.close();
}

template<typename T>
void writeFileVector(string fname, std::vector<T> vec) {
    ofstream ofs;
    ofs.open (fname, std::ofstream::out);
    for(T s : vec) {
        ofs<<s<<std::endl;
    }
    ofs.close();
}

template<typename T>
void writeFileAppendVector(string fname, std::vector<T> vec) {
    ofstream ofs;
    ofs.open (fname, std::ofstream::out | std::ofstream::app);
    for(T s : vec) {
        ofs<<s<<std::endl;
    }
    ofs.close();
}

template<typename T>
void writeFileAppendRow(string fname, std::vector<T> vec) {
    ofstream ofs;
    ofs.open (fname, std::ofstream::out | std::ofstream::app);
    for(T s : vec) {
        ofs<<s<<" ";
    }
    ofs.close();
}



namespace Utility {
using namespace SpacePrint;

// void fl();

bool isSubstring(string sub, string s){
    bool ret = false;
    if(isEmpty(sub)){
        // sub.length() == 0
        ret = true;
    }else if(sub.length() > 0 && sub.length() <= s.length()){
        // sub.length() > 0
        //
        //  abcd
        //  bc
        //   bc
        int j = 0;
        int count = 0;
        for(int i=0; i < s.length() && ret == false; i++){
            count = 0;
            j = 0;
            while(j < sub.length() && i + count < s.length() && sub[j] == s[i + count]){
                j++;
                count++;
                if(count == sub.length()){
                    ret = true;
                    break;
                }
            }
        }
    }
    return ret;
}
    
int substringFirstIndex(string sub, string s){
    int ret = -1;
    if(isEmpty(sub)){
        // sub.length() == 0
    }else if(sub.length() > 0 && sub.length() <= s.length()){
        // sub.length() > 0
        //
        //  abcd
        //  bc
        //   bc
        int j = 0;
        int count = 0;
        for(int i=0; i < s.length() && ret == -1; i++){
            count = 0;
            j = 0;
            while(j < sub.length() && i + count < s.length() && sub[j] == s[i + count]){
                j++;
                count++;
                if(count == sub.length()){
                    ret = i;
                    break;
                }
            }
        }
    }
    return ret;
}
    
// string to int
int strToInt(string s);

string fun_parent_parent();
void cut(char* pt);

int stringToInt(string s) {
    return strToInt(s);
}

int strToInt(string s) {
    return stoi(s);
}

// KEY: string to double
double stringToDouble(std::string s) {
    return stod(s);
}

// KEY string to long
long stringToLong(std::string s) {
    return stol(s);
}

// KEY string to float
float stringToFloat(std::string s) {
    return stof(s);
}

// KEY char to string
string charToString(char ch) {
    string s;
    s.push_back(ch);
    return s;
}

void begin() {
    int sz = 80;
    std::string pretty_fun = fun_parent_parent();
    int len = pretty_fun.length();
    printf ("[%s%.*s]\n", pretty_fun.c_str(), sz < len ? 0 : (int)(sz - len), PAD);
}

void begin(const char* name) {
    int sz = 80;
    printf ("[%s%.*s]\n", name, (sz < strlen(name)) ? 0 : (int)(sz - strlen(name)), PAD);
}

void end() {
    fl();
}

void cut(char* pt) {
    int k = strlen(pt) - 1;
    while(k >= 0 && pt[k] != 'v') k--;

    k >= 0 ? pt[k] = 0 : printf("Error: invalid string format.");
}

// get the name of parent function, get the name of caller function
string fun_parent() {
    void *array[10];
    size_t size;
    char **strings;

    size = backtrace (array, 10);
    strings = backtrace_symbols (array, size);

    cut(strings[1]);
    string str = std::string(strings[1] + 62);
    free(strings);
    return str;
}

// get name of caller of caller function, get the name of parent of parent function
string fun_parent_parent() {
    void *array[10];
    size_t size;
    char **strings;

    size = backtrace (array, 10);
    strings = backtrace_symbols (array, size);

    cut(strings[2]);
    string str = std::string(strings[2] + 62);
    free (strings);
    return str;
}


void print(Vector3 v) {
    v.print();
}

void print(DDLinkedList<Vector3>* ddl) {
    Node<Vector3>* curr = ddl->head;
    while(curr) {
        curr->data.print();
        curr = curr->next;
    }
}

template<typename T = int, typename U = int>
class two{
    public:
    T first;
    // T& fst = 0;
    U second;
    two(){
        // fst = 0;
    }
    two(T first, U second) {
        this -> first = first;
        this -> second = second;
    }
    public:
    string toStr(){
        string str = "(" + toStr2(first) + " " + toStr2(second) + ")";
        return str;
    }
};


}; // end class Utility




// TODO: delete it
//namespace SpaceDraw {
//class Cube {
//public:
//    float x = 1.0f;
//    float y;
//    float z;
//    float r;
//
//    GLfloat vertices[108] = { 1, 1, 1,  -1, 1, 1,  -1,-1, 1,      // v0-v1-v2 (front)
//                              -1,-1, 1,   1,-1, 1,   1, 1, 1,      // v2-v3-v0
//
//                              1, 1, 1,   1,-1, 1,   1,-1,-1,      // v0-v3-v4 (right)
//                              1,-1,-1,   1, 1,-1,   1, 1, 1,      // v4-v5-v0
//
//                              1, 1, 1,   1, 1,-1,  -1, 1,-1,      // v0-v5-v6 (top)
//                              -1, 1,-1,  -1, 1, 1,   1, 1, 1,      // v6-v1-v0
//
//                              -1, 1, 1,  -1, 1,-1,  -1,-1,-1,      // v1-v6-v7 (left)
//                              -1,-1,-1,  -1,-1, 1,  -1, 1, 1,      // v7-v2-v1
//
//                              -1,-1,-1,   1,-1,-1,   1,-1, 1,      // v7-v4-v3 (bottom)
//                              1,-1, 1,  -1,-1, 1,  -1,-1,-1,      // v3-v2-v7
//
//                              1,-1,-1,  -1,-1,-1,  -1, 1,-1,      // v4-v7-v6 (back)
//                              -1, 1,-1,   1, 1,-1,   1,-1,-1
//                            };    // v6-v5-v4
//
//    Cube() {
//    }
//    Cube(float x1, float y1, float z1, float r1=1.0) {
//        x = x1;
//        y = y1;
//        z = z1;
//        r = r1;
//    }
//
//    void draw() {
//        glEnable(GL_DEPTH_TEST);
//        glLightfv(GL_LIGHT0, GL_DIFFUSE, WHITE);
//        glLightfv(GL_LIGHT0, GL_SPECULAR, RED);
//        glMaterialfv(GL_FRONT, GL_SPECULAR, MAGENTA);
//        glMaterialf(GL_FRONT, GL_SHININESS, 10);
//        glEnable(GL_LIGHTING);
//        glEnable(GL_LIGHT0);
//
//        GLfloat lightPosition[] = {4, 3, 7, 1};
//        glLightfv(GL_LIGHT0, GL_POSITION, lightPosition);
//        glNormal3d(0, 1, 0);
//        // enble and specify pointers to vertex arrays
//        glEnableClientState(GL_VERTEX_ARRAY);
//        glVertexPointer(3, GL_FLOAT, 0, vertices);
//        glPushMatrix();
//        glDrawArrays(GL_TRIANGLES, 0, 36);
//
//        glPopMatrix();
//        glDisableClientState(GL_VERTEX_ARRAY);  // disable vertex arrays
//    }
//};  // end class Clube
//
// TODO: delete it
//class SimpleCoordinate {
//public:
//    SimpleCoordinate() {
//    }
//public:
//    void draw(float width = 1.0, int num=10) {
//        glBegin(GL_LINES);
//        float delta = width/num;
//        glColor3f(0.0f, width, 0.0f);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(-width, 0.0f, delta*i);
//            glVertex3f(width, 0.0f,  delta*i);
//        }
//
//        glColor3f(0.3f,0.7f,0.0f);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(-width, delta*i,  0.0f);
//            glVertex3f(width,  delta*i,  0.0f);
//        }
//
//        glColor3f(width, 0.0f,0.0f);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(0.0f, -width, delta*i);
//            glVertex3f(0.0f, width,  delta*i);
//        }
//
//        glColor3f(0.4f, 0.4f,0.1f);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(delta*i, -width, 0.0f);
//            glVertex3f(delta*i, width,  0.0f);
//        }
//
//        glColor3f(0.0f, 0.0f, width);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(delta*i, 0.0f, -width);
//            glVertex3f(delta*i, 0.0f, width);
//        }
//
//        glColor3f(0.0f, 0.5f, 0.5f);
//        for(int i=-num; i<=num; i++) {
//            glVertex3f(0.0f, delta*i, -width);
//            glVertex3f(0.0f, delta*i, width);
//        }
//        glEnd();
//    }
//
//}; // end class SimpleCoordinate

// TODO: delete it
//class Plane {
//public:
//    float x;
//    float y;
//    Plane(int x_ = 1.0f, int y_ = 1.0f) {
//        x = x_;
//        y = y_;
//    }
//    // draw x-y plane
//    void draw() {
//        float alpha = 0.5;
//        glBegin(GL_QUADS);
//        glColor4f(x, 0.0, 0.0, alpha);
//        glVertex3f(-x, +y, 0.0); // top left
//
//        glColor4f(0.0, y, 0.0, alpha);
//        glVertex3f(-x, -y, 0.0); // bottom left
//
//        glColor4f(0.0, 0.0, 1.0, alpha);
//        glVertex3f(+x, -y, 0.0); // bottom right
//
//        glColor4f(0.0, y, 1.0, alpha);
//        glVertex3f(+x, +y, 0.0); // top right
//        glEnd();
//    }
//}; // end class Plane
//

//TODO: delete it
//class Circle {
//    int numc;
//    double twopi;
//    double cx;
//    double cy;
//    double x;
//    double y;
//    double radius;
//    double beta = 0.0;
//public:
//    Circle(double cx, double cy) {
//        twopi = 2 * (double)M_PI;
//        this->cx = cx;
//        this->cy = cy;
//        numc = 10;
//        radius = 2;
//        beta = twopi/numc;
//    }
//
//    Circle(double cx, double cy, int numc) {
//        twopi = 2 * (double)M_PI;
//        this->numc = numc;
//        this->cx = cx;
//        this->cy = cy;
//        numc = 10;
//        radius = 2;
//        beta = twopi/this->numc;
//    }
//
//    Circle(double cx, double cy, double radius, int numc) {
//        twopi = 2 * (double)M_PI;
//        this->numc = numc;
//        this->cx = cx;
//        this->cy = cy;
//        numc = 10;
//        this->radius = radius;
//        beta = twopi/this->numc;
//    }
//    void move(double mx, double my, double mz) {
//        cx = mx;
//        cy = my;
//    }
//    void draw() {
//        GLfloat lightPosition[] = {4, 3, 7, 1};
//        glLightfv(GL_LIGHT0, GL_POSITION, lightPosition);
//        glBegin(GL_LINE_STRIP);
//        glNormal3d(0, 1, 0);
//
//        for (int i = 0; i <= numc; i++) {
//            glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE,
//                         i % 2 == 0 ? GREEN : MAGENTA);
//            x = radius*cos(beta*i);
//            y = radius*sin(beta*i);
//            glVertex3f(x + cx, y + cy, 0);
//        }
//        glEnd();
//    }
//};
//}; // end namespace PlaneSpaceDraw
//


namespace SpaceComplex {
class Complex {
public:
    double x;
    double y;

    Complex() {
        x = 0.0;
        y = 0.0;
    }
    Complex(double x_, double y_) {
        x = x_;
        y = y_;
    }

    Complex operator+(Complex& c) {
        Complex co;
        co.x = x + c.x;
        co.y = y + c.y;
        return co;
    }

    Complex operator-(Complex& c) {
        Complex co;
        co.x = x - c.x;
        co.y = y - c.y;
        return co;
    }
    Complex operator*(const Complex& c) {
        Complex c1(x*c.x - y*c.y, x*c.y + y*c.x);
        return c1;
    }

    bool operator==(const Complex& c) {
        return x == c.x && y == c.y;
    }
    Complex conjugate() {
        Complex c1(x, -y);
        return c1;
    }

    double norm(Complex& c) {
        return sqrt(c.x*c.x + c.y*c.y);
    }

    std::pair<double, double> polar() {
        double radius = norm(*this);
        double radian = acos(x/radius);
        std::pair<double, double> pair = std::make_pair(radius, radian);
        return pair;
    }


    friend Complex operator/(const Complex& c1,  const Complex& c2) {
        double xx = (c1.x*c2.x + c1.y*c2.y)/(c2.x*c2.x + c2.y*c2.y);
        double yy = (c2.x*c1.y - c1.x*c2.y)/(c2.x*c2.x + c2.y*c2.y);
        Complex c(xx, yy);
        return c;
    }


    void print() {
        printf("[%1.2f][%1.2f]\n", x, y);
    }
};

Complex rectangular(std::pair<double, double> p) {
    Complex c( p.first*cos(p.second), p.first*sin(p.second) );
    return c;
}



};

namespace MatrixVector {
class row;
class mat;

template<typename T>
T** allocateTemp(int ncol, int nrow) {
    T** arr = (T**)new T*[ncol];
    for(int c = 0; c < ncol; c++)
        arr[c] = (T*) new T[nrow];
    return arr;
}

template<typename T>
T** allocate(int ncol, int nrow) {
    T** arr = (T**)new T*[ncol];
    for(int c = 0; c < ncol; c++)
        arr[c] = (T*) new T[nrow];
    return arr;
}

template<typename T>
T** vecVecToArrArr(vector< vector<T> > vv) {
    int ncol = vv.size();
    assert(ncol > 0);
    int nrow = vv[0].size();
    T** arr = allocateTemp<T>(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            arr[c][r] = vv[c][r];
        }
    }
    return arr;
}

template<typename T>
vector< vector<T> > arrArrToVecVec(T** arr, int ncol, int nrow) {
    vector< vector<T> > v;
    for(int c=0; c<ncol; c++) {
        vector<T> v1;
        for(int r=0; r<nrow; r++) {
            v1.push_back(arr[c][r]);
        }
        v.push_back(v1);
    }
    return v;
}



mat concat(mat m1, mat m2);

class vec {
public:
    int nrow;
    int ncol;
    float **arr;
    vec();
    vec(int r, int c);
    vec(const vec &v);
    ~vec();
    vec operator+(vec v);
    bool operator==(const vec &v);
    vec& operator=(const vec &v);
    vec operator-(vec &v);
    mat operator*(row r);
    mat toMatrix();
    void print();
    row tran();
    mat multi();
    void createVec(float* arr, int len);
    vector<float> toVector();
};

class row {
public:
    int ncol;
    int nrow;
    float** arr;
public:
    row();
    ~row();
    vec tran();
    float dot(vec& r);
    row(int row);
    row(const vector<float>& vec);
    row(int col, int row);
    row(const row& r);
    void createRow(float* array, int len);
    bool operator==(const row& v);
    row operator=(const row& r);
    float operator*(vec& v);
    row operator+(row& r);
    row operator-(row& r);
    row operator/(float f);
    vector<float> toVector();

    void print();
};

vec::vec() {}

// copy constructor
vec::vec(const vec& v) {
    ncol = v.ncol;
    nrow = v.nrow;
    arr = allocate<float>(ncol, nrow);
    for(int i=0; i<ncol; i++)
        arr[i][0] = v.arr[i][0];
}

vec::vec(int col, int row=1) {
    ncol = col;
    nrow = row;
    arr = allocate<float>(ncol, nrow);
}
vec::~vec() {
    for(int i=0; i<ncol; i++) {
        if(arr[i] != NULL)
            delete[] arr[i];
    }
    if(arr != NULL)
        delete[] arr;
}

vector<float> vec::toVector() {
    vector<float> vec;
    for(int i=0; i<ncol; i++) {
        vec.push_back(arr[i][0]);
    }
    return vec;
}
vec vec::operator+(vec v) {
    assert(ncol == v.ncol);
    vec* v1 = new vec(ncol);
    for(int i=0; i<ncol; i++)
        v1->arr[i][0] = arr[i][0] + v.arr[i][0];

    return *v1;
}

vec vec::operator-(vec& v) {
    assert(ncol == v.ncol);
    vec* v1 = new vec(ncol);
    for(int i=0; i<ncol; i++)
        v1->arr[i][0] = arr[i][0] - v.arr[i][0];

    return *v1;
}


bool vec::operator==(const vec& v) {
    bool ret = true;
    for(int i=0; i<ncol && ret; i++) {
        ret = arr[i][0] == v.arr[i][0];
    }
    return ret;
}

// return reference for chaining. e.g.  (a = (b = c))
vec& vec::operator=(const vec& v) {
    assert(ncol == v.ncol);
    // same object?
    if (this != &v) {
        arr = allocate<float>(v.ncol, v.nrow);
        for(int i=0; i<ncol; i++) {
            arr[i][0] = v.arr[i][0];
        }
    }
    return *this;
}

void vec::createVec(float* array, int len) {
    for(int i=0; i<len; i++) {
        arr[i][0] = array[i];
    }
}
void vec::print() {
    for(int i=0; i<ncol; i++) {
        cout<<arr[i][0]<<std::endl;
    }
}

row::row() {
}
row::row(int col, int row) {
    ncol = col;
    nrow = row;
    arr = allocate<float>(ncol, nrow);
}
row::row(int row) {
    ncol = 1;
    nrow = row;
    arr = allocate<float>(ncol, nrow);
}
row::row(const vector<float>& vec) {
    ncol = 1;
    nrow = (int)vec.size();
    arr = allocate<float>(ncol, nrow);
    for(int i=0; i<nrow; i++) {
        arr[0][i] = vec[i];
    }
}
row::~row() {
    for(int i=0; i<ncol; i++) {
        if(arr[i] != NULL)
            delete[] arr[i];
    }
    if(arr != NULL)
        delete[] arr;
}
// copy constructor
row::row(const row& rw) {
    ncol = rw.ncol;
    nrow = rw.nrow;
    arr = allocate<float>(ncol, nrow);
    for(int r=0; r<nrow; r++) {
        arr[0][r] = rw.arr[0][r];
    }
}

row vec::tran() {
    row rw(nrow, ncol);
    for(int i=0; i<nrow; i++) {
        for(int j=0; j<ncol; j++) {
            rw.arr[i][j] = arr[j][i];
        }
    }
    return rw;
}

vec row::tran() {
    vec v = vec(nrow, ncol);
    for(int i=0; i<ncol; i++) {
        for(int j=0; j<nrow; j++) {
            v.arr[j][i] = arr[i][j];
        }
    }
    return v;
}

vector<float> row::toVector() {
    vector<float> vec;
    for(int i=0; i<nrow; i++) {
        vec.push_back(arr[0][i]);
    }
    return vec;
}
void row::print() {
    for(int i=0; i<nrow; i++) {
        cout<<arr[0][i]<<" ";
    }
    cout<<endl;
}


bool row::operator==(const row& r) {
    assert(nrow == r.nrow);
    bool ret = true;
    for(int i=0; i<nrow && ret; i++) {
        ret = arr[0][i] == r.arr[0][i];
    }
    return ret;
}

row row::operator=(const row& r) {
    ncol = r.ncol;
    nrow = r.nrow;
    arr = allocate<float>(ncol, nrow);
    for(int i=0; i<nrow; i++) {
        arr[0][i] = r.arr[0][i];
    }
    return *this;
}

void row::createRow(float* array, int len) {
    for(int i=0; i<len; i++) {
        arr[0][i] = array[i];
    }
}

class mat {
public:
    int ncol;
    int nrow;
    float** arr;
    mat();

    // copy constructor
    // called:   mat m = m1 + m2;
    // mat m(ncol, nrow);
    mat(const mat& m);
    ~mat();
    mat(const int col, const int row);
    mat(vector<vector<float>> v2);
    mat(const int col, const int row, vector<float> v);
    
    class Loader{
    public:
        mat& m;
        int ix;
        Loader(mat& m, int ix) : m(m), ix(ix) {}
        Loader operator , (float x){
            assert(ix < m.ncol * m.nrow);
            
            int u = ix / m.nrow;
            int v = ix % m.nrow;
            m.arr[u][v] = x;
            return Loader(m, ix + 1);
        }
    };
    Loader operator << (float x){
        assert(ncol > 0 && nrow > 0);
        arr[0][0] = x;
        return Loader(*this, 1);
    }
    
    vec getVec(int n);
    row getRow(int n);
    vec operator*(vec& v);
    mat operator*(mat& m);
    mat operator=(const mat& m);
    bool operator==(const mat& m);
    mat operator+(mat m);
    mat operator-(const mat& m);
    mat operator/(float f);
    mat operator*(float f);
    mat removeRow(int index);
    mat removeVec(int index);
    mat insertVecNext(int index, vec v);
    mat insertVecPrevious(int index, vec v);
    mat clone();
    mat concat(mat m);


    mat rowMultiScala(int index, float f);
    mat vecMultiScala(int index, float f);
    mat swapRow(int inx1, int inx2);

    // TODO: move to as normal functions
    void geneMat(int init);
    
    // TODO: the method should be removed, there is geneMatRandom(int ncol, nrow, int fst, int snd)
    void geneMatRandom(int fst, int snd);
    
    void print();
    string toStr();
    void print2();
    void zero();
    void identity();
    // check if the matrix is an identity matrix
    bool isId();
    // check if the matrix is zero matrix
    bool isZero();
    float det();
    
    mat id();

    mat subMatrix(int colIndex, int rowIndex);

    mat transpose();

    /*
    1 2 3 block(1, 1, 1, 2) => 5 6
    4 5 6
    */
    mat block(int colIndex, int clen, int rowIndex, int rlen);
    mat take(int len);
    mat drop(int len);
    mat init();
    mat tail();
//        mat addVecAfter(int index, float f);
//        mat rowAdd(int index, float f);
//        mat rowSub(int index, float f);
//        mat vecSub(int index, float f);
//        mat subVec(int index, vec v);
//        mat subRow(int index, row r);
//        mat upperTri();
//        mat lowerTri();
//        bool isSymmetric();
//        float trace();
//        mat subMatrix(int col, int row);

};

mat::mat() {
    ncol = 0;
    nrow = 0;
    arr = NULL;
    cout<<"empty matrix";
}
void mat::zero() {
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            arr[c][r] = 0;
        }
    }
}

mat mat::id() {
    assert(ncol == nrow);
    mat m(nrow, ncol);
    m.identity();
    return m;
}

// identity matrix
void mat::identity() {
    assert(ncol == nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            if( c == r)
                arr[c][r] = 1;
            else
                arr[c][r] = 0;
        }
    }
}

// m1 = [[1, 2],
//       [3, 4]]
// m2 = [[1, 0],
//       [0, 1]]
// m = concat(m1, m2)
// m = [[1, 2, 1, 0],
//      [3, 4, 0, 1]]
mat::mat(const mat& m) {
    ncol = m.ncol;
    nrow = m.nrow;
    arr = allocate<float>(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            arr[c][r] = m.arr[c][r];
        }
    }
}
mat::mat(int col, int row) {
    ncol = col;
    nrow = row;
    arr = allocate<float>(ncol, nrow);
}
mat::mat(vector<vector<float>> v2){
    assert(v2.size() > 0 && v2[0].size() > 0);
    ncol = v2.size();
    nrow = v2[0].size();
    arr = allocate<float>(ncol, nrow);
    for(int c=0; c<ncol; c++){
        for(int r=0; r<nrow; r++)
            arr[c][r] = v2[c][r];
    }
}
mat::mat(int col, int row, vector<float> v){
    ncol = col; 
    nrow = row;
    assert(ncol >= 0 && nrow >= 0 && v.size() == ncol*nrow);
    arr = allocate<float>(ncol, nrow);
    for(int c=0; c<ncol; c++){
        for(int r=0; r<nrow; r++)
            arr[c][r] = v[nrow*c + r];
    }
}

mat::~mat() {
    for(int i=0; i<ncol; i++)
        delete[] arr[i];

    delete[] arr;
}

mat vec::toMatrix() {
    mat m(ncol, 1);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m.arr[c][r] = arr[c][r];
        }
    }
    return m;
}

void mat::print() {
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            cout<<arr[c][r]<<" ";
        }
        cout<<std::endl;
    }
}

string mat::toStr() {
    string ret = "";
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            ret += toStr2(arr[c][r]) + " ";
            // cout<<arr[c][r]<<" ";
        }
        ret += "\n";
    }
    return ret;
}

void mat::print2() {
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            printf("[%*f]", 10, arr[c][r]);
        }
        printf("\n");
    }
}

vec mat::getVec(int n) {
    assert(0 <= n && n < nrow);
    vec v(ncol);
    for(int i=0; i<ncol; i++) {
        v.arr[i][0] = arr[i][n];
    }
    return v;
}

row mat::getRow(int index) {
    row r(1, nrow);
    for(int i=0; i<nrow; i++) {
        r.arr[0][i] = arr[index][i];
    }
    return r;
}

mat mat::removeVec(int index) {
    assert(nrow > 0 && index >= 0 && index < nrow);
    mat m(ncol, nrow - 1);
    for(int c=0; c<ncol; c++) {
        int k = 0;
        for(int r=0; r<nrow; r++) {
            if(index != r) {
                m.arr[c][k] = arr[c][r];
                k++;
            }
        }
    }
    return m;
}

mat mat::insertVecNext(int index, vec v) {
    assert(ncol == v.ncol && 0 <= index && index < nrow);
    mat mv = v.toMatrix();
    mat left = take(index + 1);
    mat right = drop(index + 1);
    mat cat = left.concat(mv).concat(right);
    return cat;
}

/*
Add vector to the left of index
Shoud I change the name to insertVecLeft()
'left' is shorter than 'previous'
*/
mat mat::insertVecPrevious(int index, vec v) {
    assert(ncol == v.ncol && 0 <= index && index < nrow);
    mat mv = v.toMatrix();
    if (index == 0) {
        return mv.concat(*this);
    } else {
        mat left = take(index);
        mat right = drop(index);
        return left.concat(mv).concat(right);
    }
}

mat mat::removeRow(int index) {
    assert(ncol > 0 && index >= 0 && index < ncol);
    mat m(ncol - 1, nrow);
    int k = 0;
    for(int c=0; c<ncol; c++) {
        if(index != c) {
            for(int r=0; r<nrow; r++) {
                m.arr[k][r] = arr[c][r];
            }
            k++;
        }
    }
    return m;
}
mat mat::operator+(mat m) {
    mat m1(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m1.arr[c][r] = arr[c][r] + m.arr[c][r];
        }
    }
    return m1;
}

mat mat::operator-(const mat &m) {
    mat m1(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m1.arr[c][r] = arr[c][r] - m.arr[c][r];
        }
    }
    return m1;
}
mat mat::operator/(float f) {
    mat m1(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m1.arr[c][r] = arr[c][r]/f;
        }
    }
    return m1;
}

mat mat::operator*(float f) {
    mat m1(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m1.arr[c][r] = arr[c][r]*f;
        }
    }
    return m1;
}


mat mat::subMatrix(int cStart, int rStart) {
    assert(cStart < ncol && rStart < nrow);
    int ncolt = ncol - cStart;
    int nrowt = nrow - rStart;
    mat m(ncolt, nrowt);
    for(int c=0; c<m.ncol; c++) {
        for(int r=0; r<m.nrow; r++) {
            m.arr[c][r] = arr[c + cStart][r + rStart];
        }
    }
    return m;
}

mat mat::take(int len) {
    mat m = block(0, ncol, 0, len);
    return m;
}

mat mat::drop(int len) {
    mat m = block(0, ncol, len, nrow - len);
    return m;
}
mat mat::init() {
    assert(ncol > 0 && nrow > 0);
    mat m = take(nrow - 1);
    return m;
}
mat mat::tail() {
    assert(ncol > 0 && nrow > 0);
    mat m = drop(1);
    return m;
}


mat mat::block(int cIndex, int clen, int rIndex, int rlen) {
    mat m(clen, rlen);
    for(int c=0; c<m.ncol; c++) {
        for(int r=0; r<m.nrow; r++) {
            m.arr[c][r] = arr[cIndex + c][rIndex + r];
        }
    }
    return m;
}

mat mat::transpose() {
    mat m(nrow, ncol);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m.arr[r][c] = arr[c][r];
        }
    }
    return m;
}

// two empty matrices are equal
bool mat::operator==(const mat& m) {
    bool ret = ncol == m.ncol && nrow == m.nrow ? true : false;
    for(int c=0; c<ncol && ret; c++) {
        for(int r=0; r<nrow && ret; r++) {
            ret = arr[c][r] == m.arr[c][r];
        }
    }
    return ret;
}

void mat::geneMat(int init) {
    int n = init;
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            arr[c][r] = n;
            n++;
        }
    }
}

void mat::geneMatRandom(int fst, int snd){
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            arr[c][r] = random(snd - fst) + fst;
        }
    }
}


mat vec::operator*(row rw) {
    // M(m,n) M(n,h) => M(m,h)
    assert(nrow == rw.ncol);
    mat m(ncol, rw.nrow);
    for(int c=0; c<m.ncol; c++) {
        for(int r=0; r<m.nrow; r++) {
            m.arr[c][r] = arr[c][0]*rw.arr[0][r];
        }
    }
    return m;
}

mat vec::multi() {
    mat m(2, 2);
    return m;
}

vec mat::operator*(vec &v) {
    assert(nrow == v.ncol);
    vec v1(ncol);
    for(int c=0; c<ncol; c++) {
        v1.arr[c][0] = getRow(c)*v;
    }
    return v1;
}
mat mat::operator*(mat &m) {
    assert(nrow == m.ncol);
    mat m1(ncol, m.nrow);
    m1.zero();
    for(int r=0; r<nrow; r++) {
        m1 = m1 + getVec(r) * m.getRow(r);
    }
    return m1;
}
//  c = (a = b) ?
mat mat::operator=(const mat& m) {
    if (this != &m) {
        ncol = m.ncol;
        nrow = m.nrow;
        arr = allocate<float>(ncol, nrow);
        for(int c=0; c<ncol; c++) {
            for(int r=0; r<nrow; r++) {
                arr[c][r] = m.arr[c][r];
            }
        }
    }
    return *this;
}

mat mat::clone() {
    mat m(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m.arr[c][r] = arr[c][r];
        }
    }
    return m;
}
// concate two matrices
mat mat::concat(mat m) {
    assert(ncol == m.ncol);
    mat m1(ncol, nrow + m.nrow);
    for(int c=0; c < ncol; c++) {
        for(int r=0; r < nrow + m.nrow; r++) {
            if (r < nrow)
                m1.arr[c][r] = arr[c][r];
            else {
                m1.arr[c][r] = m.arr[c][r - nrow];
            }
        }
    }
    return m1;

//        nrow = nrow + m.nrow;
//        arr = allocate<float>(ncol, nrow);
//        for(int c=0; c < ncol; c++){
//            for(int r=0; r < nrow; r++){
//                arr[c][r] = m1.arr[c][r];
//            }
//        }
}

mat mat::rowMultiScala(int index, float f) {
    mat m(ncol, nrow);
    for(int c=0; c < ncol; c++) {
        for(int r=0; r < nrow; r++) {
            if(index != c) {
                m.arr[c][r] = arr[c][r];
            } else {
                m.arr[c][r] = arr[c][r]*f;
            }
        }
    }
    return m;
}

mat mat::vecMultiScala(int index, float f) {
    mat m(ncol, nrow);
    for(int c=0; c < ncol; c++) {
        for(int r=0; r < nrow; r++) {
            if(r != index) {
                m.arr[c][r] = arr[c][r];
            } else {
                m.arr[c][r] = arr[c][r]*f;
            }
        }
    }
    return m;
}

mat mat::swapRow(int inx1, int inx2){
    assert(0 <= inx1 && inx1 < ncol && 0 <= inx2 && inx2 < ncol);
    row r1 = getRow(inx1);
    row r2 = getRow(inx2);
    mat m(ncol, nrow);
    for(int c=0; c < ncol; c++) {
        for(int r=0; r < nrow; r++) {
            if(c == inx1){
                m.arr[c][r] = r2.arr[0][r];
            }else if(c == inx2){
                m.arr[c][r] = r1.arr[0][r];
            }else{
                m.arr[c][r] = arr[c][r];
            }
        }
    }
    return m;
}

float row::dot(vec& v) {
    assert(nrow == v.ncol);
    float ret = 0.0;
    for(int i=0; i<nrow; i++) {
        ret += arr[0][i]*v.arr[i][0];
    }
    return ret;
}



float row::operator*(vec& v) {
    return dot(v);
}
row row::operator+(row& r) {
    assert(nrow == r.nrow);
    row rw(ncol, nrow);

    for(int i=0; i<nrow; i++) {
        rw.arr[0][i] = arr[0][i] + r.arr[0][i];
    }
    return rw;
}
row row::operator-(row& r) {
    assert(nrow == r.nrow);
    row rw(ncol, nrow);

    for(int i=0; i<nrow; i++) {
        rw.arr[0][i] = arr[0][i] - r.arr[0][i];
    }
    return rw;
}
row row::operator/(float f) {
    row rw(ncol, nrow);
    for(int i=0; i<nrow; i++) {
        rw.arr[0][i] = arr[0][i]/f;
    }
    return rw;
}
mat concat(mat m1, mat m2) {
    assert(m1.ncol == m2.ncol);

    mat m(m1.ncol, m1.nrow + m2.nrow);
    for(int c=0; c < m1.ncol; c++) {
        for(int r=0; r < m1.nrow + m2.nrow; r++) {
            if (r < m1.nrow)
                m.arr[c][r] = m1.arr[c][r];
            else {
                m.arr[c][r] = m2.arr[c][r - m1.nrow];
            }
        }
    }
    return m;
}

mat identity(int n) {
    mat m(n, n);
    for(int c=0; c<n; c++) {
        for(int r=0; r<n; r++) {
            if( c == r)
                m.arr[c][r] = 1;
            else
                m.arr[c][r] = 0;
        }
    }
    return m;
}

vec backwardSubstitute(mat a, vec b) {
    assert(a.ncol == a.nrow && a.nrow == b.ncol);

    vec x(a.ncol);
    for(int c=a.ncol - 1; c >= 0; c--) {
        float s = 0.0;
        for(int r= a.nrow - 1; r >=c; r--) {
            if(r == c) {
                x.arr[r][0] = (b.arr[c][0] - s)/a.arr[c][r];
            } else {
                s += a.arr[c][r]*x.arr[r][0];
            }
        }
    }
    return x;
}

mat ltri(mat m, int index) {
    assert(0 <= index && index < m.nrow);
    mat im = identity(m.nrow);
    vec v = m.getVec(index);

    // assume m.arr[index][index] != 0
    for(int c=m.ncol-1; c > index; c--) {
        im.arr[c][index] = - m.arr[c][index]/m.arr[index][index];
    }
    return im;
}

std::pair<mat, mat> utri(mat m) {
    mat m4;
    mat m2 = m4;
    mat m1 = m;
    mat l = identity(m.ncol);
    for(int r = 0; r < m.nrow - 1; r++) {
        mat mt = ltri(m1, r);
        pl("what");
        l = mt*l;  // L_k..L2L1
        m1 = l*m;  // L_k..L2 (L1 A)
    }
    return std::make_pair(l, l*m);
}

mat geneMat(int ncol, int nrow, int init){
    mat m(ncol, nrow);
    int count = init;
    for(int c=0; c<ncol; c++){
        for(int r=0; r<nrow; r++){
            m.arr[c][r] = count;
            count++;
        }
    }
    return m;
}
    
mat geneMatRandom(int ncol, int nrow, int fst, int snd){
    mat m(ncol, nrow);
    for(int c=0; c<ncol; c++) {
        for(int r=0; r<nrow; r++) {
            m.arr[c][r] = random(snd - fst) + fst;
        }
    }
    return m;
}
    

    
    
}; // end MatrixVector

namespace SpaceVector4 {
class Vector4 {
    float column[4];
public:
    const float e1[4] = {1.0f, 0.0f, 0.0f, 0.0f};
    const float e2[4] = {0.0f, 1.0f, 0.0f, 0.0f};
    const float e3[4] = {0.0f, 0.0f, 1.0f, 0.0f};
    const float e4[4] = {0.0f, 0.0f, 0.0f, 1.0f};
public:

    Vector4() {
        column[0] = column[1] = column[2] = column[3] = 0.0f;
    }

    Vector4(const Vector4& other) {
        this->column[0] = other.column[0];
        this->column[1] = other.column[1];
        this->column[2] = other.column[2];
        this->column[3] = other.column[3];
    }

    Vector4(float x, float y, float z, float w = 1.0f) {
        this->column[0] = x;
        this->column[1] = y;
        this->column[2] = z;
        this->column[3] = w;
    }
    Vector4(const float arr[4]) {
        this->column[0] = arr[0];
        this->column[1] = arr[1];
        this->column[2] = arr[2];
        this->column[3] = arr[3];
    }
    Vector4(float arr[4]) {
        this->column[0] = arr[0];
        this->column[1] = arr[1];
        this->column[2] = arr[2];
        this->column[3] = arr[3];
    }
    bool operator==(const Vector4& rhs) {
        bool b0 = column[0] == rhs.column[0];
        bool b1 = column[1] == rhs.column[1];
        bool b2 = column[2] == rhs.column[2];
        bool b3 = column[3] == rhs.column[3];
        return (b0 && b1 && b2 && b3);
    }
    Vector4& operator=(const Vector4& rhs) {
        this->column[0] = rhs.column[0];
        this->column[1] = rhs.column[1];
        this->column[2] = rhs.column[2];
        this->column[3] = rhs.column[3];
        return *this;
    }

    Vector4 operator+(Vector4& rhs) {
        Vector4 v;
        v.column[0] = this->column[0] + rhs.column[0];
        v.column[1] = this->column[1] + rhs.column[1];
        v.column[2] = this->column[2] + rhs.column[2];
        v.column[3] = 0.0f;
        return v;
    }

    Vector4 operator-(Vector4& rhs) {
        Vector4 v;
        v.column[0] = this->column[0] - rhs.column[0];
        v.column[1] = this->column[1] - rhs.column[1];
        v.column[2] = this->column[2] - rhs.column[2];
        v.column[3] = 0.0f;
        return v;
    }

    Vector4 operator/(float n) {
        Vector4 v;
        v.column[0] = column[0]/n;
        v.column[1] = column[1]/n;
        v.column[2] = column[2]/n;
        return v;
    }

    float dot(Vector4& rhs) {
        Vector4 v;
        v.column[0] = column[0] * rhs.column[0];
        v.column[1] = column[1] * rhs.column[1];
        v.column[2] = column[2] * rhs.column[2];
        v.column[3] = column[3] * rhs.column[3];
        return v[0] + v[1] + v[2] + v[3];
    }

    float cross(Vector4& rhs) {
        return 0.0;
    }

    Vector4 normal() {
        Vector4 v;
        float norm = column[0]*column[0] +
                     column[1]*column[1] +
                     column[2]*column[2];
        float n = sqrtf(norm);
        return (*this)/n;
    }

    // [] can't modify member variables
    const float& operator[](int index) const {
        return column[index];
    }

    // [] can modify member variables
    float& operator[](int index) {
        return column[index];
    }

    void pp() {
        print();
    }
    void print() {
        printf("x=[%1.2f]\ny=[%1.2f]\nz=[%1.2f]\nw=[%1.2f]\n\n", column[0], column[1], column[2], column[3]);
    }
};
};



// C++14 only
// http://localhost:8080/snippet?id=s+c%2B%2B
#if __cplusplus >= 201402

std::pair<string, string> split(string s, int inx) {
    unsigned long len = s.length();
    if(0 <= inx && inx <= len) {
        string prefix = s.substr(0, inx);
        string suffix = s.substr(inx, len);
        return make_pair(prefix, suffix);
    } else if(inx < 0) {
        return make_pair("", s);
    } else {
        return make_pair(s, "");
    }
}

template<typename T, typename Fun>
static vector<vector<T>> splitWhen(Fun f, const vector<T>& vec) {
    vector<vector<T>> ret;
    vector<T> v2;
    for(auto const& e : vec) {
        if(!f(e)) {
            v2.push_back(e);
        }else{
            if(len(v2) > 0){
                vector<T> tmp(v2);
                ret.push_back(tmp);
                v2.clear();
            }
        }
    }
    if(len(v2) > 0){
        vector<T> tmp(v2);
        ret.push_back(tmp);
    }
    return ret;
}




string removeIndex(string s, int inx) {
    std::pair<string, string> p = split(s, inx);
    string prefix = p.first;
    string suffix = p.second;
    string tail = suffix.substr(1, suffix.length());
    return prefix + tail;
}
template<typename T>
vector<T> removeIndex(vector<T>& vec, int inx){
    vector<T> ret(vec); // clone vec
    ret.erase(ret.begin() + inx);
    return ret;
}

template<typename T>
vector<T> removeIndexRange(vector<T>& vec, int fromInx, int toInx){
    vector<T> ret(vec); // clone vec
    // Error if vector is removed from the smaller index
    for(int i=toInx; i >= fromInx; i--){
        ret.erase(ret.begin() + i);
    }
    return ret;
}




namespace SpaceMatrix4 {
using namespace SpaceVector4;
using namespace Utility;

class Matrix4 {
    Vector4 mat[4];
public:
    Matrix4() {
        mat[0] = mat[1] = mat[2] = mat[3] = Vector4();
    }

    // copy constructor
    Matrix4(const Matrix4& matrix) {
        mat[0] = matrix[0];
        mat[1] = matrix[1];
        mat[2] = matrix[2];
        mat[3] = matrix[3];
    }

    Matrix4(Vector4 v0, Vector4 v1, Vector4 v2, Vector4 v3) {
        mat[0] = v0;
        mat[1] = v1;
        mat[2] = v2;
        mat[3] = v3;
    }
    Matrix4(float m[16]) {
        Vector4 v1({ m[0],   m[1],   m[2],   m[3]});
        Vector4 v2({ m[4],   m[5],   m[6],   m[7]});
        Vector4 v3({ m[8],   m[9],   m[10],  m[11]});
        Vector4 v4({ m[12],  m[13],  m[14],  m[15]});
        mat[0] = v1;
        mat[1] = v2;
        mat[2] = v3;
        mat[3] = v4;
    }
    Matrix4& operator=(const Matrix4& matrix) {
        mat[0] = matrix[0];
        mat[1] = matrix[1];
        mat[2] = matrix[2];
        mat[3] = matrix[3];
        return *this;
    }

    bool operator==(const Matrix4& matrix) {
        bool b0 = mat[0] == matrix[0];
        bool b1 = mat[1] == matrix[1];
        bool b2 = mat[2] == matrix[2];
        bool b3 = mat[3] == matrix[3];
        return (b0 && b1 && b2 && b3);
    }
    // overload [] , const => member variables CAN NOT be modified
    const Vector4& operator[](int index) const {
        return mat[index];
    }

    // overload [] , no const => member variables CAN be modified
    Vector4& operator[](int index) {
        return mat[index];
    }
    Matrix4 operator+(Matrix4& rhs) {
        Matrix4 m;
        m[0] = mat[0] + rhs[0];
        m[1] = mat[1] + rhs[1];
        m[2] = mat[2] + rhs[2];
        m[3] = mat[3] + rhs[3];
        return m;
    }
    Matrix4 operator-(Matrix4& rhs) {
        Matrix4 m;
        m[0] = mat[0] - rhs[0];
        m[1] = mat[1] - rhs[1];
        m[2] = mat[2] - rhs[2];
        m[3] = mat[3] - rhs[3];
        return m;
    }

    // column vector
    Vector4 operator*(Vector4& vect4) {
        Vector4 row0(mat[0][0], mat[1][0], mat[2][0], mat[3][0]);
        Vector4 row1(mat[0][1], mat[1][1], mat[2][1], mat[3][1]);
        Vector4 row2(mat[0][2], mat[1][2], mat[2][2], mat[3][2]);
        Vector4 row3(mat[0][3], mat[1][3], mat[2][3], mat[3][3]);

        row0.print();
        row1.print();
        row2.print();
        row3.print();

        vect4.print();
        cout<<"["<<row0.dot(vect4)<<"]"<<std::endl;
        cout<<"["<<row1.dot(vect4)<<"]"<<std::endl;
        cout<<"["<<row2.dot(vect4)<<"]"<<std::endl;
        cout<<"["<<row3.dot(vect4)<<"]"<<std::endl;

        Vector4 v(row0.dot(vect4), row1.dot(vect4), row2.dot(vect4), row3.dot(vect4));
        v.print();
        return v;
    }

    // column vector
    Matrix4 operator*(Matrix4& matrix) {
        Matrix4 m;
        m[0] = (*this)*matrix[0];
        m[1] = (*this)*matrix[1];
        m[2] = (*this)*matrix[2];
        m[3] = (*this)*matrix[3];
        return m;
    }
    // [1, 2, 3, 0] => point at infinite
    // [1, 2, 3, 1] => normal point
    Matrix4 translate(float x, float y, float z) {
        Matrix4 m;
        mat[0][0] = 1;
        mat[1][1] = 1;
        mat[2][2] = 1;
        mat[3][3] = 1;

        mat[3][0] = x;
        mat[3][1] = y;
        mat[3][2] = z;
        mat[3][3] = 1;
        return *this;
    }
    Matrix4 identity() {
        mat[0][0] = 1;
        mat[1][1] = 1;
        mat[2][2] = 1;
        mat[3][3] = 1;
        return *this;
    }

    void pp() {
        print();
    }
    void print() {
        printf("[%1.2f][%1.2f][%1.2f][%1.2f]\n", mat[0][0], mat[1][0], mat[2][0], mat[3][0]);
        printf("[%1.2f][%1.2f][%1.2f][%1.2f]\n", mat[0][1], mat[1][1], mat[2][1], mat[3][1]);
        printf("[%1.2f][%1.2f][%1.2f][%1.2f]\n", mat[0][2], mat[1][2], mat[2][2], mat[3][2]);
        printf("[%1.2f][%1.2f][%1.2f][%1.2f]\n", mat[0][3], mat[1][3], mat[2][3], mat[3][3]);
        fl();
    }
};
};

namespace Lambda {


template<typename T>
T head(vector<T> vec) {
    assert(vec.size() > 0);
    return vec[0];

}

template<typename T>
T last(vector<T> vec) {
    assert(vec.size() > 0);
    return vec.back();
}


template<typename T>
static vector<T> reverse(vector<T> vec) {
    vector<T> v;
    for(auto it=vec.rbegin(); it != vec.rend(); it++)
        v.push_back(*it);

    return v;
}

template<typename Fun, typename T>
static vector<T> takeWhile(Fun f, vector<T> vec) {
    vector<T> ret;
    for(auto const& n : vec) {
        if(f(n))
            ret.push_back(n);
        else
            break;
    }
    return ret;
}

template<typename Fun, typename T>
static vector<T> dropWhile(Fun f, vector<T> vec) {
    vector<T> ret;
    bool flag = false;
    for(auto const& n : vec) {
        if(!f(n) || flag) {
            ret.push_back(n);
            flag = true;
        }
    }
    return ret;
}

template<typename T>
static vector<T> tail(vector<T> vec) {
    vector<T> v;
    for(int i=1; i<vec.size(); i++) {
        v.push_back(vec.at(i));
    }
    return v;
}

template<typename Fun, typename T>
static T foldl(Fun f, T acc, vector<T> vec) {
    vector<T> vec1;
    T tmpAcc = acc;
    for(auto ite = vec.begin(), end = vec.end(); end != ite; ite++) {
        tmpAcc = f(tmpAcc, *ite);
    }
    return tmpAcc;
}


template<typename Fun, typename T>
static T foldr(Fun f, T acc, const vector<T>& vec) {
    vector<T> vec1;
    T tmpAcc = acc;
    for(int i=vec.size() - 1; i >= 0; i--) {
        tmpAcc = f(vec[i], tmpAcc);
    }
    return tmpAcc;
}

template<class T>
static vector<T> init(const vector<T>& vec) {
    assert(vec.size() > 0);
    vector<T> v;
    for(int i=0; i < vec.size() - 1; i++) {
        v.push_back(vec[i]);
    }
    return v;
}

template<typename Fun, typename T>
static vector<T> map(Fun f, vector<T> vec) {
    vector<T> vec1;
    for(auto const& v : vec) {
        vec1.push_back(f(v));
    }
    return vec1;
}

int min(int a, int b) {
    return a < b ? a : b;
}

template<typename Fun, typename T>
static vector<T> zipWith(Fun f, vector<T> v1, vector<T> v2) {
    vector<T> vec;
    int len1 = unsignToInt(v1.size());
    int len2 = unsignToInt(v2.size());
    for(int i=0; i<min(len1, len2); i++) {
        for(int j=0; j<min(len1, len2); j++) {
            vec.push_back(f(v1[i], v2[j]));
        }
    }
    return vec;
}

template<typename T>
static vector<T> flat(const vector<vector<T>>& vec2){
    vector<T> vec;
    for(auto& v : vec2){
        vec = vec + v;
    }
    return vec;
}

template<typename T, typename Fun>
static vector<T> filter(Fun f, const vector<T>& vec) {
    vector<T> vec1;
    for(auto const& v : vec) {
        if(f(v)) {
            vec1.push_back(v);
        }
    }
    return vec1;
}

// take n char from a string
static string take(int n, string s) {
    return s.substr(0, n);
}

template<class T>
static vector<T> take(int n, const vector<T>& vec) {
    vector<T> vec1;
    int c = 0;
    for(auto const& v : vec) {
        if (c < n) {
            vec1.push_back(v);
        } else {
            break;
        }
        c++;
    }
    return vec1;
}

template<class T>
static vector<T> drop(int n, vector<T> vec) {
    vector<T> vec1;
    int c = 0;
    for(auto const& v : vec) {
        if (c >= n)
            vec1.push_back(v);
        c++;
    }
    return vec1;
}

template<typename Fun, typename T>
static vector<T> mergeSortListLam(Fun f, vector<T> v1, vector<T> v2) {
    vector<T> empty;
    if (v1.size() == 0)
        return v2;
    if (v2.size() == 0)
        return v1;

    vector<T> vh1 = take(1, v1);
    vector<T> vh2 = take(1, v2);
    if (f(vh1.at(0), vh2.at(0))) {
        return vh1 + mergeSortListLam(f, tail(v1), v2);
    } else {
        return vh2 + mergeSortListLam(f, v1, tail(v2));
    }
}

template<class T>
static vector<T> mergeSortList(vector<T> v1, vector<T> v2) {
    vector<T> empty;
    if (v1.size() == 0)
        return v2;
    if (v2.size() == 0)
        return v1;

    vector<T> vh1 = take(1, v1);
    vector<T> vh2 = take(1, v2);
    if (vh1.at(0) < vh2.at(0)) {
        return vh1 + mergeSortList(tail(v1), v2);
    } else {
        return vh2 + mergeSortList(v1, tail(v2));
    }
}

template<class T>
static vector<T> mergeSort(vector<T> v1) {
    if(v1.size() > 1) {
        vector<T> left = take(v1.size()/2, v1);
        vector<T> right= drop(v1.size()/2, v1);
        vector<T> subl = mergeSort(left);
        vector<T> subr = mergeSort(right);
        return mergeSortList(subl, subr);
    } else {
        return v1;
    }
}

}; // end namespace Lambda

template<class T>
vector<T> interleave(const vector<T>& v1, const vector<T>& v2) {
    int len1 = v1.size();
    int len2 = v2.size();
    vector<T> vec;
    for(int i=0; i<std::min(len1, len2); i++) {
        vec.push_back(v1[i]);
        vec.push_back(v2[i]);
    }
    return vec;
}

string takeFileName(string path) {
    string s;
    // it works for MacOS so far, not windows
    vector<string> vec = splitStr(path, "/");
    if(vec.size() > 0) {
        s = Lambda::last(vec);
    }
    return s;
}

string takeDirectory(string path) {
    vector<string> vec = splitStr(path, "/");
    vector<string> vecDir = Lambda::init(vec);
    vector<string> pathvec;
    int len = vecDir.size();
    vector<string> vv = interleave(vecDir, repeatVec(len, c2s("/")));
    if(!(Lambda::head(vec).empty() && vecDir.size() == 1)) {
        pathvec = Lambda::init(vv);
    } else {
        pathvec = vv;
    }

    return Lambda::foldr([](auto x, auto y) {
        return x + y;
    }, c2s(""), pathvec);
}

string takeExtension(string path) {
    auto vec = splitStr(path, ".");
    string ret;
    if(vec.size() > 1) {
        string dot(".");
        ret = dot + Lambda::last(vec);
    }
    return ret;
}

string join(vector<string>& vec) {
    return trimRight(Lambda::foldr([](auto& x, auto&y) {
        return x + " " + y;
    }, c2s(""), vec));
}

string joinNewLine(vector<string> vec) {
    return Lambda::foldr([](auto& x, auto&y) {
        return x + " " + y;
    }, c2s("\n"), vec);
}

vector<string> readFile(string fname) {
    std::ifstream file(fname);
    vector<string> retVec;
    if (file.is_open()) {
        std::string line;
        while(getline(file, line)) {
            retVec.push_back(line);
        }
        file.close();
    }
    return retVec;
}

void listDirToFile(char * path, char* toFile, vector<string>& vec) {
    DIR * d = opendir(path);
    if( d != NULL) {
        struct dirent* dir;
        while ((dir = readdir(d)) != NULL) {
            if(dir->d_type != DT_DIR) {
                char fpath[1000];
                sprintf(fpath, "%s/%s", path, dir->d_name);
                vec.push_back(c2s(fpath));
                if(vec.size() > 1000) {
                    writeFileAppendVector(c2s(toFile), vec);
                    vec.clear();
                }
            } else if(dir -> d_type == DT_DIR && strcmp(dir->d_name,".")!=0 && strcmp(dir->d_name,"..")!=0 ) {
                char fpath[1000];
                sprintf(fpath, "%s/%s", path, dir->d_name);
                listDirToFile(fpath, toFile, vec);
            }
        }
        closedir(d);
    }
}


void listDirToVec(char * path, char* toFile, vector<string>& vec) {
    listDirToFile(path, toFile, vec);
}



vector<string> recurveDirBoost(char* dir) {
    vector<string> retVec;
    fs::path p(dir);
    bool b = exists(p);

    if(!exists(p) || !is_directory(p)) {
        std::cout << p << " is not a path\n";
        return retVec;
    }

    fs::recursive_directory_iterator begin(p), end;
    std::vector<fs::directory_entry> v(begin, end); // vector created here
    for(auto& f: v) {
        retVec.push_back(f.path().string());
    }

    return retVec;
}




void recurveDirToFile(char * path, char* toFile, vector<string>& vec) {
    DIR * d = opendir(path);
    if( d != NULL) {
        struct dirent* dir;
        while ((dir = readdir(d)) != NULL) {
            if(dir-> d_type != DT_DIR) {
                char fpath[500];
                sprintf(fpath, "%s/%s", path, dir->d_name);
                vec.push_back(c2s(fpath));
                if(vec.size() > 1000) {
                    writeFileAppendVector(c2s(toFile), vec);
                    vec.clear();
                }
            } else if(dir -> d_type == DT_DIR && strcmp(dir->d_name,".")!=0 && strcmp(dir->d_name,"..")!=0 ) {
                char fpath[500];
                sprintf(fpath, "%s/%s", path, dir->d_name);
                listDirToFile(fpath, toFile, vec);
            } else {
                vec.push_back("Unknown_type");
            }
        }
        if(vec.size() > 0) {
            writeFileAppendVector(c2s(toFile), vec);
            vec.clear();
        }
        closedir(d);
    }
}

void listDirAll(char * path) {
    DIR * d = opendir(path);
    if(d != NULL) {
        struct dirent* dir;
        while ((dir = readdir(d)) != NULL) {
            if(dir-> d_type != DT_DIR) { // if the type is not directory just print it with blue
                char dpath[1000]; // here I am using sprintf which is safer than strcat
                sprintf(dpath, "%s/%s", path, dir->d_name);
                printf("%s\n", dpath);
            } else if(dir -> d_type == DT_DIR && strcmp(dir->d_name,".")!=0 && strcmp(dir->d_name,"..")!=0 ) {
                printf("%s\n", dir->d_name); // print its name in green
                char dpath[1000]; // here I am using sprintf which is safer than strcat
                sprintf(dpath, "%s/%s", path, dir->d_name);
                printf("%s\n", dpath);
                listDirAll(dpath); // recall with the new path
            } else {
                pl("Unknown_type");
            }
        }
        closedir(d); // finally close the directory
    }
}

void listDirToVector(char * path, vector<string>& vec) {
    DIR * d = opendir(path);
    if( d != NULL) {
        struct dirent* dir;
        while ((dir = readdir(d)) != NULL) {
            if(dir-> d_type != DT_DIR) {
                char fullpath[1000];
                sprintf(fullpath, "%s/%s", path, dir->d_name);
                vec.push_back(c2s(fullpath));
            } else if(dir -> d_type == DT_DIR && strcmp(dir->d_name,".")!=0 && strcmp(dir->d_name,"..")!=0 ) {
                char fullpath[1000];
                sprintf(fullpath, "%s/%s", path, dir->d_name);
                listDirToVector(fullpath, vec);
            } else {
                vec.push_back("Unknown_type");
            }
        }
        closedir(d);
    }
}


string getPWD() {
    string ret;
    char cwd[PATH_MAX];
    if (getcwd(cwd, sizeof(cwd)) != NULL) {
        ret = cStringToString(cwd);
    } else {
        perror("getcwd() error");
    }
    return ret;
}

string intToStringLen(int num, int maxlen) {
    char buffer[255];
    sprintf(buffer, "%d", num);
    string zero = Lambda::take(maxlen - strlen(buffer), repeat(num, "0"));
    return zero + c2s(buffer);
}

std::pair<vector<string>, vector<string>> parseLambda(string s) {
    vector<string> vecStr = splitStrRegex(s, "\\s+");
    vector<string> vec = Lambda::filter([](auto& x) {
        return trim(x).length() > 0;
    }, vecStr);

    vector<string> varv;
    auto f = [](auto x) {
        return x == c2s("->");
    };
    std::pair<vector<string>, vector<string>> pair = breakVec(f, vec);
    return std::make_pair(pair.first, Lambda::tail(pair.second));
}




namespace Algorithm {
using namespace Lambda;

void merge(int* arr, int lo, int mid, int hi) {
    int len = hi - lo + 1;
    int* list = new int[len];
    int i = lo;
    int j = mid + 1;
    int k = 0;
    while(i < mid + 1 || j < hi + 1) {
        if (i >= mid + 1) {
            list[k] = arr[j];
            j++;
        } else if(j >= hi + 1) {
            list[k] = arr[i];
            i++;
        } else if (arr[i] < arr[j]) {
            list[k] = arr[i];
            i++;
        } else {
            list[k] = arr[j];
            j++;
        }
        k++;
    }
    for(int i=0; i<len; i++) {
        arr[lo + i] = list[i];
    }
    delete[] list;
}

void mergeSort(int* arr, int lo, int hi) {
    if(lo < hi) {
        int m = (lo + hi)/2;
        mergeSort(arr, lo, m);
        mergeSort(arr, m+1, hi);
        merge(arr, lo, m, hi);
    }
}

// TODO: add test cases
void swap(int array[], int i, int j) {
    int tmp = array[i];
    array[i] = array[j];
    array[j] = tmp;
}

    
int partition(int array[], int lo, int hi) {
    if(hi > lo) {
        int p = array[hi];
        int big = lo;
        for(int i=lo; i<=hi; i++) {
            if(array[i] <= p)     {
                swap(array, i, big);
                if(i < hi)
                    big++;
            }
        }
        return big;
    }
    return -1;
}

// TODO: add test cases
void quickSortArr(int array[], int lo, int hi) {
    if(hi > lo) {
        int pivot = partition(array, lo, hi);
        quickSortArr(array, lo, pivot-1);
        quickSortArr(array, pivot+1, hi);
    }
}

// TODO: add lambda fun for more general sorting
template<typename T>
static vector<T> quickSort(vector<T> vec) {
    if (vec.size() <= 1)
        return vec;
    else {
        T pivot = head(vec);
        vector<T> rest = drop(1, vec);
        vector<T> left =  filter<T>([pivot](auto x) {
            return x < pivot;
        }, rest);
        vector<T> right = filter<T>([pivot](auto x) {
            return x >= pivot;
        }, rest);
        return quickSort(left) + con(pivot, quickSort(right));
    }
}

template<typename Fun, typename T>
static vector<T> mergeSortList(Fun f, vector<T> v1, vector<T> v2) {
    vector<T> empty;
    if (v1.size() == 0)
        return v2;
    if (v2.size() == 0)
        return v1;

    vector<T> vh1 = take(1, v1);
    vector<T> vh2 = take(1, v2);
    if (f(vh1.at(0), vh2.at(0))) {
        return vh1 + mergeSortList(f, tail(v1), v2);
    } else {
        return vh2 + mergeSortList(f, v1, tail(v2));
    }
}

}
#endif

//string toStr(int e){
//    string ret;
//    return ret;
//}
//string toStr(float e){
//    string ret;
//    return ret;
//}
//string toStr(double e){
//    string ret;
//    return ret;
//}
//string toStr(long e){
//    string ret;
//    return ret;
//}

using namespace MatrixVector;
void writeFile(string fname, mat& m) {
    for(int i=0; i<m.nrow; i++) {
        vector<float> vec = m.getVec(i).toVector();
        writeFileAppendVector(fname, vec);
    }
}

//void writeFileRow(string fname, row& r){
//    auto f = [](auto& x) { return toStr(s);};
//    vector<float> vec = {0.1, 1.2};
//    vector<string> s1 = Lambda::map( f, vec);
//}


#endif