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 <sstream>
#include <string>
#include <sstream>
#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 "/Users/aaa/myfile/bitbucket/clib/AronCLibNew.h"
 
// #include "/Users/cat/myfile/bitbucket/cpplib/AronPrint.h"
// #include <experimental/optional>
// #include "/Library/Developer/CommandLineTools/usr/include/c++/v1/experimental/optional"
// #include "/usr/local/include/boost/operations.hpp"
 
#include <vector>
 
// KEY: math library
// URL: https://openframeworks.cc/documentation/glm/
#include "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};
#include <boost/filesystem.hpp>
 
#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";
 
// KEY: namespace forward declaration
/*
namespace MatrixVector {
  class row;
  class mat;
};
*/
 
namespace fs = boost::filesystem;
using namespace std;
// using namespace AronGeometry;
 
int len(string s);
 
template<typename T>
int len(vector<T> v);
 
template<typename T>
string toStr(const T& s);
 
template<typename T>
string toStr2(T& s);
 
template<typename T>
string toStrVar(T& s);
 
template<typename T>
vector<vector<T>> transpose(vector<vector<T>> mat);
 
// NOTE: forward template class NO default type
// ERROR: template<typename T = int , typename U = int> 
template<typename T, typename U>
class two;
 
template<typename T>
bool compareVector(vector<T> v1, vector<T> v2);
 
template<typename T>
bool compareVector(vector<vector<T>> v1, vector<vector<T>> v2);
 
template<typename T>
vector<vector<T>> arrToVec2(T* arr, int len, int width);
 
void printArr2(float* arr, int len, int width);
 
 
 
int unsignToInt(size_t u);
 
string padRight(string s, string pad, int maxLen);
string padLeft(string s, string pad, int maxLen);
 
template<typename T>
two<T, T> divMod(T a, T b);
 
void newline();
void nl();
int max(int a, int b, int c);
int max(int a, int b, int c, int d);
int max(vector<int> v);
 
string reverseStr(string s);
 
int maxConsecutiveOne(vector<int> vec);
int strToInt(string s);
int numFromBase(string fromBase, string input);
 
vector<int> range(int a, int b);
vector<int> range(int a, int b, int stride);
 
template<typename T>
vector<T> reverseVec(vector<T> vec);
 
std::map<char, int> hexMap();
 
string baseToBase(string fromBase, string toBase, string input);
 
inline uint64_t nowMilliseconds();
inline uint64_t nowMicroseconds();
inline uint64_t nowNanoseconds();
 
template<typename T>
vector<T> repeatVec(int n, T t);
template<typename T>
vector<T> concatVec(vector<T> v1, vector<T> v2);
// vector<T> concatVec(vector<T>& v1, vector<T>& v2);
 
vector<int> randomNanoSec(int n, int max);
 
template<typename T>
T dotProduct(vector<T> v1, vector<T> v2);
 
template<typename T>
T det2(vector<vector<T>> v2);
 
template<typename T>
vector<T> subtract(vector<T> v1, vector<T> v2);
 
 
void testfun()
{
}
 
void testfun1() {  
}
 
void testfun2 ( ) {  
}
 
 
inline uint64_t nowMilliseconds() {
    return std::chrono::duration_cast<std::chrono::milliseconds>
              (std::chrono::high_resolution_clock::now().time_since_epoch()).count();
}
 
inline uint64_t nowMicroseconds() {
    return std::chrono::duration_cast<std::chrono::microseconds>
              (std::chrono::high_resolution_clock::now().time_since_epoch()).count();
}
inline uint64_t nowNanoseconds() {
    return std::chrono::duration_cast<std::chrono::nanoseconds>
              (std::chrono::high_resolution_clock::now().time_since_epoch()).count();
}  
 
 
// string to int
int strToInt(string s) {
    return stoi(s);
}
 
// KEY: get size of types
void getTypeSize(){
    cout<<"int           => " <<sizeof(int)<<endl;
    cout<<"long          => " <<sizeof(long)<<endl;
    cout<<"float         => " <<sizeof(float)<<endl;
    cout<<"double        => " <<sizeof(double)<<endl;
    cout<<"size_t        => " <<sizeof(size_t)<<endl;
    cout<<"unsigned int  => " <<sizeof(unsigned int)<<endl;
    cout<<"unsigned long => " <<sizeof(unsigned long)<<endl;
    cout<<"long long     => " <<sizeof(long long)<<endl;
}
 
int maxConsecutiveOne(vector<int> vec){
  int max = 0;
  int prev = -1;
  int curr = -1;
  for(int i = 0; i < len(vec); i++){
    if (vec[i] == 1){
      if(prev == -1){
        // Previous position
        prev = i;
        curr = i;
      }else{
        // Current position
        curr = i;
      }
    }else{
      if(prev != -1 && curr != -1){
        int dist = curr - prev + 1;
        if(dist > max){
          max = dist;
        }
      }
      // Reset positions
      prev = -1;
      curr = -1;
    }
  }
  
  if(prev != -1 && curr != -1){
    int dist = curr - prev + 1;
    if(dist > max){
      max = dist;
    }
  }
  return max;
}
 
int maxConsecutiveChar(char* arr, int size, char c){
  int max = 0;
  int prev = -1;
  int curr = -1;
  for(int i = 0; i < size; i++){
    if(arr[i] == c){
      if(prev == -1){
        prev = i;
        curr = i;
      }else{
        curr = i;
      }
    }else{
      if(prev != -1){
        int dist = curr - prev + 1;
        max = dist > max ? dist : max;
        prev = -1;
        curr = -1;
      }
    }
  }
  if(prev != -1){
    int dist = curr - prev + 1;
    max = dist > max ? dist : max;
  }
  return max;
}
 
std::map<char, int> hexMap(){
  char c1 = '0';
  char c2 = 'a';
  int k = 0;
  
  std::map<char, int> map;
  for(int i = 0; i < 16; i++){
    if(i < 10){
      map[c1 + i] = i;
    }
    else{
      map[c2 + k] = i;
      k++;
    }
  }
  return map;
}
 
std::map<int, char> hexMapRev(){
  std::map<int, char> revMap;
  std::map<char, int> map = hexMap();
  std::map<char, int>::iterator it = map.begin();
  while(it != map.end()){
    char key = it -> first;
    int value = it -> second;
    revMap[value] = key;
    it++;
  }
  return revMap;
}
 
int numFromBase(string fromBase, string input){
  int k = 0;
  int sum = 0;
  std::map<char, int> map = hexMap();
  for(auto i : reverseVec(range(0, input.size()))){
    int ii = input.size() - 1 - i;
    if('0' <= input[i] && input[i] <= '9' || 'a' <= input[i] && input[i] <= 'z'){
      int value = map[input[i]];
      if(value != NULL){
        sum += value * pow(strToInt(fromBase), ii);
      }
    }
  }
  return sum;
}
 
string baseToBase(string fromBase, string toBase, string input){
  string ret = "";
  int sum = numFromBase(fromBase, input);
  std::map<int, char> map = hexMapRev();
  int toBaseNum = strToInt(toBase);
  if(sum == 0){
    ret += '0';
  }else{
    while(sum > 0){
      int q = sum / toBaseNum;
      int r = sum % toBaseNum;
      ret = map[r] + ret;
      sum = q;
    }      
  }
  return ret;
}
 
vector<int> toBinary(int n){
  vector<int> vec;
  string s = baseToBase("10", "2", toStr2(n));
  for(int i = 0; i < len(s); i++){
    if (s[i] == '0'){
      vec.push_back(0);
    }else if(s[i] == '1'){
      vec.push_back(1);
    }else{
      printf("ERROR: Invalid Char");
      exit(1);
    }
  }
  return vec;
}
 
template<typename T>
vector<vector<T>> arrToVec2(T* arr, int len, int width){
    vector<vector<T>> vv;
 
    int q = len / width;
    int r = len % width;
 
    for(int i = 0; i < q; i++){
        vector<T> v;
        for(int j = 0; j < width; j++){
            v.push_back(arr[width * i + j]);
        }
        vv.push_back(v);
    }
    if(0 < r){
        vector<T> w;
        for(int i = 0; i < width; i++){
            if(i < r){
                w.push_back(arr[width * q + i]);
            }else{
                w.push_back(0.0);
            }
        }
        vv.push_back(w);
    }
    return vv;
}
 
 
template<typename T>
two<T, T> divMod(T a, T b){
  two<T, T> t(a / b, a % b);
  return t;
}
 
 
int max(int a, int b, int c){
    return max(max(a, b), c);
}
 
int max(int a, int b, int c, int d){
    return max(max(max(a, b), c), d);
}
 
int max(vector<int> v){
  if(v.size() > 0){
    int m = v[0];
    for(auto x : v){
      m = max(m, x);
    }
    return m;
  }else{
    printf("ERROR: vector can not be empty\n");
    std::exit(EXIT_FAILURE);
  }
}
 
int abs_int(int a){
  std::abs(a);
}
float abs_f(float a){
  std::abs(a);
}
long abs_l(long a){
  std::abs(a);
}
long long abs_ll(long long a){
  std::abs(a);
}
double abs_d(double a){
  std::abs(a);
}
unsigned int abs_ui(unsigned int a){
  return 0;
}
 
string reverseStr(string s){
    std::string rs = s;
    std::reverse(rs.begin(), rs.end());
    return rs;
}
 
 
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();
}
 
void newline(){
  printf("\n");
}
 
 
 
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 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();
}
 
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, "");
    }
}
 
vector<int> sumVec(vector<int> v1, vector<int> v2){
  vector<int> vec;
  int diff = v2.size() > v1.size() ? v2.size() - v1.size() : v1.size() - v2.size();
  vector<int> mv = v2.size() > v1.size() ? v2 : v1;
  vector<int> iv = v2.size() > v1.size() ? v1 : v2;
  vector<int> vz = repeatVec(diff, 0);
  vector<int> w1 = mv;
  vector<int> w2 = concatVec(vz, iv);
  vector<int> ux = repeatVec(len(mv) + 1, 0);
  
  int c = 0;
  for(int i = len(w1) - 1; i >= 0; i--){
    int s = w1[i] + w2[i] + c;
    ux[i + 1] = s % 10;
    c = s / 10;
  }
  ux[0] = c;
 
  return ux;
}
 
 
 
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;
}
 
 
string repeat(int n, string s);
string cStringToString(const char* pt);
 
template<typename T>
void pl(T t);
 
void print(glm::vec3 v3);
void print(glm::vec2 v2);
string removeAdjacentDuplicateChar(string s, char c);
string concatValidPath(vector<string> vec);
 
 
 
template<typename T>
void print(std::pair<T, T> pair);
 
template<typename T>
void print(vector<T> vec);
 
template<typename T>
void printVec(vector<T> vec);
 
 
template<typename T>
void print(vector< vector<T> > vec);
 
template<typename T>
void print(T t);
 
template<typename T>
void print(string msg, T t);
 
void nl();
 
template<typename T>
void printLnNb(T t);
 
template<typename T>
vector<T> operator+(vector<T> v1, vector<T> v2);
 
string operator+(const char* chars, string s);
 
string operator+(string s, const char* chars);
 
template<typename T>
vector<T> con(T t, vector<T> v);
 
namespace AronPrint {
using namespace std;
 
void ppf(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) {
  std::cout<<"call me maybe"<<endl;
  std::cout<<t1<<t2<<endl;
}
 
template<typename T>
void pp(T t1, T t2, T t3) {
  std::cout<<t1<<t2<<t3<<endl;
}
void ppf(const char* format, ...) {
  std::cout<<"printf format"<<endl;
  printf(format);
}
 
template<typename T, typename U>
void pp(std::pair<T, U> p){ 
  pp(toStr(p));
}
 
void pp(int n) {
  std::cout<<"["<<n<<"]"<<endl;
}
 
void pp(float n) {
  std::cout<<"["<<n<<"]"<<endl;
}
 
void pp(long n) {
  std::cout<<"["<<n<<"]"<<endl;
}
void pp(unsigned long n) {
  std::cout<<"["<<n<<"]"<<endl;
}
 
void pp(string msg0, string msg1){
  std::cout<<msg0<<msg1<<std::endl;
}
void pp(string msg0, float num){
  std::cout<<msg0<<num<<std::endl;
}
void pp(string msg0, int num){
  std::cout<<msg0<<num<<std::endl;
}
void pp(string msg0, long num){
  std::cout<<msg0<<num<<std::endl;
}
void pp(string msg0, double num){
  std::cout<<msg0<<num<<std::endl;
}
void pp(string msg0, char* charPt){
  std::cout<<msg0<<charPt<<std::endl;
}
 
/*
void fl() {
  newline();
  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);
  cout<<endl;
  printLnNb(repeat(pad, "-") + s + repeat(pad, "-"));
}
 
} // end AronPrint
 
 
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 printTick() {
    endTime = clock();
    clock_t diff = endTime - startTime;
    std::cout<<"Tick=["<<diff<<"] CPU ticks"<<std::endl;
  }
};
 
 
template<typename T = int, typename U = int>
class two{
 public:
  T first;
  T x;
  // T& fst = 0;
  U second;
  U y;
  two(){
    // fst = 0;
  }
  two(T first, U second) {
    this -> x = this -> first = first;
    this -> y = this -> second = second;
  }
 public:
  string toStr(){
    return "(" + toStr2(x) + "," + toStr2(y) + ")";
  }
  string toString(){
    return toStr();
  }
};
 
 
 
namespace AronLambda {
template<typename T>
T head(vector<T> vec) {
    assert(vec.size() > 0);
    return vec[0];
 
}
 
char head(string s){
  if(len(s) > 0)
    return s[0];
  else
    throw std::runtime_error("string is empty.");
}
 
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 T>
static vector<T> rev(vector<T> vec) {
  return reverse(vec);
}
 
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 Fun>
static string takeWhile(Fun f, string str) {
  string ret = "";
  int inx = -1;
  for(int i = 0; i < len(str); i++){
    if(!f(str[i])){
      inx = i;
      break;
    }
  }
  if(inx != -1)
    ret = str.substr(0, inx);
  
  return ret;
}
 
template<typename Fun>
static string dropWhile(Fun f, string str) {
  string ret = "";
  int inx = -1;
  for(int i = 0; i < len(str); i++){
    if(!f(str[i])){
      inx = i;
      break;
    }
  }
  if(inx != -1)
    ret = str.substr(inx, len(str));
  
  return ret;
}
 
template<typename Fun>
static two<string, string> splitWhileStr(Fun f, string str){
  two<string, string> tw;
  
  int ln = len(str);
  int i = 0;
  for(i = 0; i < ln; i++){
    if(!f(str[i])){
      break;
    }
  }
  string s1 = str.substr(0, i);
  string s2 = str.substr(i, ln);
  tw.x = s1;
  tw.y = s2;
  
  return tw;
}
 
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) {
    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) {
    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 A>
static vector<A> map(Fun f, vector<A> vec) {
    vector<A> vec1;
    for(auto& v : vec) {
        vec1.push_back(f(v));
    }
    return vec1;
}
 
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(vector<vector<T>>& vec2){
    vector<T> vec;
    for(vector<T>& 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);
}
 
/*
string cStringToString(const char* pt) {
    std::string s(pt);
    return s;
}
*/
 
static string take(int n, const char* str){
  string s = cStringToString(str);
  return take(n, s);
}
 
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 AronLambda
 
template<typename T, typename U>
class two;
 
using namespace AronLambda;
 
 
 
// using namespace MySpace;
 
// namespace fs = boost::filesystem;
 
// using namespace Algorithm;
// using namespace AronClass;
 
 
 
 
template<typename T>
void pl(T t) {
    print(t);
    newline();
}
 
// TODO: Rename a better name
template<typename T>
string toStr2(T& s) {
    std::ostringstream out;
    out<<s;
    return out.str();
}
 
template<typename T>
string toStrVar(T& s) {
    std::ostringstream out;
    out<<s;
    return out.str();
}
 
 
template<typename T>
void swap(vector<T>& vec, int i, int j);
 
 
int random(int n);
string cStringToString(const char* pt);
 
void nl(){
    printf("\n");
}
 
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;
}
 
 
long lcgRandom(long x0){
  long a = 1664525;
  long c = 1013904223;
  long m = (long)pow((double)2, (double)32);
  long x1 = (a * x0 + c) % m;
  return x1;
}
 
 
vector<int> randomIntVec(int count, int max){
    vector<int> v;
    long x0 = 22337737;
    for(int i = 0; i < count; i++){
        int x1 = lcgRandom(x0) % INT_MAX;
        v.push_back(x1 % max);
        x0 = (long)x1;
    }
    return v;
}
 
vector<int> randomNanoSec(int n, int max){
    vector<int> v;
    for(int i = 0; i < n; i++){
        auto n = nowNanoseconds();
        v.push_back(n % max);
    }
    return v;
}
 
 
 
int random(int n) {
  long x0 = 399;
  return lcgRandom(x0) % 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;
}
 
vector<int>* nPrime(int n){
  vector<int>* vec = new vector<int>();
  if(n > 0){
    vec -> push_back(2);
    int c = 1;
    int num = 3;
    while(c < n){
      bool isPrime = true;
      for(auto const& p : *vec){
        if (num % p == 0){
          isPrime = false;
          break;
        }
      }
      if(isPrime){
        vec -> push_back(num);
        c++;
      }
      num++;
    }
  }
  return std::move(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);
}
 
const char* strToCharArray(string s){
  return s.c_str();
}
 
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 str2(T& s) {
    return toStr2(s);
}
 
string trimLeft(string str) {
    string s = str;
    s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](char& x) {
        return !std::isspace(x);
    }));
    return s;
}
 
string trimRight(string str) {
    string s = str;
    s.erase(std::find_if(s.rbegin(), s.rend(), [](char& 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 charPtrToString(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 = [](double xx, int nn, double cc) {
        double nb = intToDouble(nn);
        return pow(xx, nb) - cc;
    };
 
    // The first derivative of f(x)
    auto df = [](double xx, int nn, double 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;
}
 
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;
}
 
template<typename T>
void add(vector<T>& vec, T t){
  vec.push_back(t);
}
 
int compareCString(char* s1, char* s2) {
    return strcmp(s1, s2);
}
 
int compareString(string s1, string s2) {
    return s1.compare(s2);
}
 
template<typename T>
bool compareVector(vector<T> v1, vector<T> v2){
    if(v1.size() == v2.size()){
        for(int i = 0; i < v1.size(); i++){
            if(v1[i] != v2[i]){
                return false;
            }
        }
        return true;
    }
    return false;
}
 
template<typename T>
bool compareVector(vector<vector<T>> v1, vector<vector<T>> v2){
    if(v1.size() == v2.size()){
        for(int i = 0; i < v1.size(); i++){
          if(!compareVector(v1[i], v2[i])){
            return false;
          }
        }
        return true;
    }
    return false;
}
 
 
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;
}
 
int inx(int cIx, int rIx, int width){
    return cIx * width + rIx;
}
 
void tranmat4(float arr[16]){
    int width = 4;
    for(int c = 0; c < width; c++){
        for(int r = c; r < width; r++){
            float tmp = arr[inx(c, r, width)];
            arr[inx(c, r, width)] = arr[inx(r, c, width)];
            arr[inx(r, c, width)] = tmp;
        }
    }
}
 
 
// 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(two<T, T> t) {
  cout<<"("<<t.x<<" "<<t.y<<")";
}
 
 
template<typename T>
void print(vector<T> vec) {
    for(auto const& v: vec) {
        std::cout<<v<<" ";
    }
    
    std::cout<<std::endl;
}
 
template<typename T>
void printVec(vector<T> vec) {
    print(vec);
}
 
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<<"];"
  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 ppL(T t) {
    printLn(t);
}
 
 
void printArr(float arr[10], int n, int nCol){
  for(int i = 0; i < n; i++){
    if((i + 1) % nCol != 0){
      cout<<arr[i]<<" ";
    }else{
      cout<<arr[i]<<" "<<endl;
    }
  }
}
 
template<typename T>
vector<T> concatVec(vector<T> v1, vector<T> v2){
// vector<T> concatVec(vector<T>& v1, vector<T>& v2){
  return v1 + v2;
}
 
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());
}
 
 
 
 
template<typename T>
int length(vector<T> v){
    return unsignToInt(v.size());
}
 
int len(string s){
    return length(s);
}
 
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;
}
 
 
template<typename T>
vector<T> reverseVec(vector<T> vec){
  vector<T> v(vec);
  std::reverse(v.begin(), v.end());
  return v;
}
 
 
// 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;
}
 
vector<int> range(int a, int b, int stride) {
    vector<int> vec;
    for(int i=a; i<=b; i += stride)
        vec.push_back(i);
    return vec;
}
 
// generate vector double
vector<int> geneVector(int a, int b) {
    int len = b - a;
    std::vector<int> vec;
    for(int i=0; i<=len; i++) {
        vec.push_back(a + i);
    }
    return vec;
}
 
bool isBST(Node<int>* curr, Node<int>* &prev){
  if(curr != NULL){
    if(!isBST(curr -> left, prev))
      return false;
    if(prev != NULL && prev -> data > curr -> data)
      return false;
    prev = curr;
    if(!isBST(curr -> right, curr))
      return false;
  }
  return true;
}
 
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>
int partitionT(T array[], int lo, int hi) {
    if(hi > lo) {
        T p = array[hi];
        int bInx = lo;
        for(int i=lo; i<=hi; i++) {
            if(array[i] <= p)     {
                if(i != bInx){
                    swapT(array, i, bInx);
                }
 
                if(i < hi){
                    bInx++;
                }
            }
        }
        return bInx;
    }
    return lo;
}
 
// 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>
int partitionInline(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){
                if (i != bigInx){
                    swap(vec, i, bigInx);
                }
                if(bigInx < hi){
                  bigInx++;
                }
            }
        }
    }
    return bigInx;
}
 
template<typename T>
void quickSortVec(vector<T>& vec, int lo, int hi){
    if(lo < hi){
        int pInx = partitionInline(vec, lo, hi);
        quickSortVec(vec, lo, pInx - 1);   
        quickSortVec(vec, pInx + 1, hi);
    }
}
 
template<typename T>
void quickSortT(T array[], int lo, int hi) {
    if(hi > lo) {
        int pInx = partitionT(array, lo, hi);
        quickSortT(array, lo, pInx - 1);
        quickSortT(array, pInx + 1, hi);
    }
}
 
std::string replace(std::string str, regex rex, std::string rep){
    return std::regex_replace (str, rex, rep);
}
 
bool isEmptyLine(string s){
  regex rx("\\s+");
  string str = replace(s, rx, "");
  return len(str) == 0;
}
 
// get environment variable
string getEnv(string env){
  const char* pt = std::getenv(s2c(env));
  return c2s(pt);
}
 
template<typename T>
vector<T> getCol(int ix, vector<vector<T>> v){
  vector<vector<T>> tv = transpose(v);
  if(0 <= ix && ix < tv.size()){
    return tv[ix];
  }else{
    printf("ERROR: Invalid index=%d\n", ix);
  }
}
 
 
string concatValidPath(vector<string> vec){
  string s = AronLambda::foldl([](string x, string y) { return x + "/" + y;}, c2s(""), vec);
  return removeAdjacentDuplicateChar(s, '/');
}
 
string removeAdjacentDuplicateChar(string s, char c){
  string ret;
  if(len(s) > 0) {
    ret.push_back(s[0]);
    for (int i = 0; i + 1 < len(s); i++) {
      if (s[i] != s[i + 1]) {
        ret.push_back(s[i + 1]);
      }else{
        if(c != s[i + 1])
          ret.push_back(s[i + 1]);
      }
    }
  }
  return ret;
}
 
#if 1 
namespace SpaceTest {
using namespace AronPrint;
 
    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("True=> " + 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>
    static void t(vector<T> v1, vector<T> v2, string s = "") {
        if(compareVector(v1, v2)){
            pp("True => " + s);
        }else{
            printVec(v1);
            printVec(v2);
            pp("False => " + s);
        }
        nl();
    }
  
    template<typename T>
    static void t(vector<vector<T>> v1, vector<vector<T>> v2, string s = "") {
        if(compareVector(v1, v2)){
            pp("True => " + s);
        }else{
            print(v1);
            print(v2);
            pp("False => " + s);
        }
        nl();
    }
};  // end SpaceTest 
#endif
 
// all the geometry functions here
 
namespace AronGeometry {
using namespace AronPrint;
 
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;
    }
 
 
    // 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>
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);
    }
};
 
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;
}
 
} // send AronGeometry
 
 
 
 
 
 
namespace Utility {
using namespace AronPrint;
 
// 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 fun_parent_parent();
void cut(char* pt);
 
  
int stringToInt(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;
}
 
// KEY: padding string
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);
}
 
// KEY: padding string
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;
    }
}
 
};      // 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
 
 
 
 
 
 
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;
}
 
 
 
}; // end SpaceComplex
 
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(const vector<float> &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();
    vector<float> toVec();
    int size();
};
 
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();
    vector<float> toVec();
 
    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(const vector<float> &v){
  ncol = v.size();
  nrow = 1;
  arr = allocate<float>(ncol, nrow);
  cout<<"nrow="<<nrow<<endl;
  cout<<"ncol="<<ncol<<endl;
  for(int i=0; i < ncol; i++){
    arr[i][0] = v[i];
  }
}
  
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;
    }
}
 
int vec::size(){
    return ncol;
}
vector<float> vec::toVector() {
    vector<float> vec;
    for(int i=0; i<ncol; i++) {
        vec.push_back(arr[i][0]);
    }
    return vec;
}
vector<float> vec::toVec() {
  return toVector();
}
  
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 = ncol == v.ncol;
  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;
}
  
vector<float> row::toVec() {
  return toVector();
}
  
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];
    }
}
 
// Thu 15 Jun 10:35:28 2023  
// TODO: Use template<typename T>
//       T** arr
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);
    }
 
    vector<vector<float>> toVecVec();
    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(){
  vector<int> vmax;
  for(int i = 0; i < nrow; i++){
    vector<int> vm;
    for(int j = 0; j < ncol; j++){
      vm.push_back(len(toStrVar(arr[j][i])));
    }
    vmax.push_back(max(vm));
  }
 
  vector<vector<string>> vv;
  string pad = " ";
  for(int i = 0; i < nrow; i++){
    vector<string> vs;
    for(int j = 0; j < ncol; j++){
      string s = padRight(toStrVar(arr[j][i]), pad, vmax[i] + 1);
      vs.push_back(s);
    }
    vv.push_back(vs);
  }
  vector<vector<string>> rv = ::transpose(vv);
  ::print(rv);
}
 
  
 
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;
}
  vector<vector<float>> mat::toVecVec(){
    vector<vector<float>> vv;
    for(int c=0; c<ncol; c++) {
      vector<float> tv;
        for(int r=0; r<nrow; r++) {
          tv.push_back(arr[c][r]);
        }
        vv.push_back(tv);
    }
    return vv;
  }
  
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){
    long r0 = 9947;
    long max = (long)(snd - fst);
    for(int c=0; c<ncol; c++) {
        for(int i=0; i<nrow; i++) {
          long r1 = lcgRandom(r0) + fst;
          arr[c][i] = r1;
          r0 = r1;
        }
    }
}
 
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;
            // m.arr[c][r] = 9.7; 
        }
    }
    m.print();
    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]);
    }
};
}; // end SpaceVector4
 
// C++14 only
// http://localhost:8080/snippet?id=s+c%2B%2B
#if __cplusplus >= 201402
 
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();
    }
};
};  // end SpaceMatrix4
 
template<class T>
vector<T> interleave(const vector<T>& v1, const vector<T>& v2) {
  vector<T> vec;
  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;
}
 
using namespace AronPrint;
 
string takeFileName(string path) {
  string s;
  // it works for MacOS so far, not windows
  vector<string> vec = splitStr(path, "/");
  if(vec.size() > 0) {
    s = AronLambda::last(vec);
  }
  return s;
}
 
string takeDirectory(string path) {
    vector<string> vec = splitStr(path, "/");
    vector<string> vecDir = AronLambda::init(vec);
    vector<string> pathvec;
    int len = vecDir.size();
    vector<string> vv = interleave(vecDir, repeatVec(len, c2s("/")));
    if(!( AronLambda::head(vec).empty() && vecDir.size() == 1)) {
        pathvec = AronLambda::init(vv);
    } else {
        pathvec = vv;
    }
 
    return AronLambda::foldr([](auto x, auto y) {
        return x + y;
    }, c2s(""), pathvec);
}
 
 
 
template<typename T>
vector<T> interleave(vector<T> v1, vector<T> v2, vector<T> v3){
  vector<T> vec;
  int vlen = std::min(std::min(len(v1), len(v2)), len(v3));
  for(int i = 0; i < vlen; i++){
    vec.push_back(v1[i]);
    vec.push_back(v2[i]);
    vec.push_back(v3[i]);
  }
  return vec;
}
 
 
string takeExtension(string path) {
    auto vec = splitStr(path, ".");
    string ret;
    if(vec.size() > 1) {
        string dot(".");
        ret = dot + AronLambda::last(vec);
    }
    return ret;
}
 
string join(vector<string>& vec) {
    return trimRight( AronLambda::foldr([](auto x, auto y) {
        return x + " " + y;
    }, c2s(""), vec));
}
 
string joinNewLine(vector<string> vec) {
    return AronLambda::foldr([](auto x, auto y) {
        return x + " " + y;
    }, c2s("\n"), vec);
}
 
string concatR(vector<string>& vec, string s) {
    return AronLambda::foldr([&s](auto x, auto y) {
        return x + s + y;
    }, c2s(""), vec);
}
 
string concatL(vector<string>& vec, string s) {
    return AronLambda::foldl([&s](auto x, auto y) {
        return x + s + y;
    }, c2s(""), 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;
}
 
 
vector<string> readFileBackward(string fname, int nline){
    vector<string> vec;
    char c;
    int maxChar = 1000;
    std::ifstream myFile(fname,std::ios::ate);
    std::streampos size = myFile.tellg();
    char* pt = (char*)malloc(maxChar*sizeof(char));
    int nx = 0;
    int k = 0;
    for(int i=1;i<=size && k < nline;i++){
        myFile.seekg(-i,std::ios::end);
        myFile.get(c);
        if(c != '\n'){
            pt[nx] = c; 
            nx++;
            if(nx >= maxChar){
                printf("ERROR: nx > maxChar\n");
                printf("ERROR: readFileBackward(string fname, int nline)\n");
                break;
            }
        }else{
            pt[nx] = '\0'; 
            vec.push_back(reverseStr(c2s(pt)));
            free(pt);
            pt = (char*)malloc(maxChar*sizeof(char));
            nx = 0;
            k++;
 
            // free(pt) ?
        }
    }
    return reverse(vec);
}
 
 
string readFileStream(const char* fpath){
  std::string vertexCode;
  std::ifstream vShaderFile;
  // ensure ifstream objects can throw exceptions:
  vShaderFile.exceptions (std::ifstream::failbit | std::ifstream::badbit);
  try {
    // open files
    vShaderFile.open(fpath);
    std::stringstream vShaderStream;
    // read file's buffer contents into streams
    vShaderStream << vShaderFile.rdbuf();
    // close file handlers
    vShaderFile.close();
    // convert stream into string
    vertexCode   = vShaderStream.str();
  }
  catch (std::ifstream::failure& e) {
    std::cout << "ERROR::SHADER::FILE_NOT_SUCCESFULLY_READ" << std::endl;
  }
  return vertexCode;
}
 
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(fs::directory_entry& 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 = AronLambda::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 = AronLambda::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, AronLambda::tail(pair.second));
}
 
#endif
 
namespace Algorithm {
 
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));
    }
}
 
}  // end Algorithm
 
//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);
    }
}
 
template<typename T>
T dotProduct(vector<T> v1, vector<T> v2){
  T n;
  for(int i = 0; i < min(len(v1), len(v2)); i++){
    n += v1[i] * v2[i];
  }
  return n;
}
 
template<typename T>
T det2(vector<vector<T>> v2){
  if (len(v2) == 2 && len(v2[0]) == 2){
    T a = v2[0][0];
    T b = v2[0][1];
    T c = v2[1][0];
    T d = v2[1][1];
    return a*d - b*c;
  }else{
    pp("Invalid Input.");
  }
}
 
template<class T>
T det2(vector<T> v1, vector<T> v2){
  if (len(v1) == 2 && len(v2) == 2){
    T a = v1[0];
    T b = v1[1];
    T c = v2[0];
    T d = v2[1];
    return a*d - b*c;
  }else{
    pp("Invalid Input.");
  }
}
 
 
template<class T>
vector<T> crossProduct(vector<T> v1, vector<T> v2){
  if(len(v1) == len(v2) && len(v1) == 3){
    T x;
    {
      vector<T> v11 = removeIndex(v1, 0);
      vector<T> v22 = removeIndex(v2, 0);
      x = det2<T>(v11, v22);
    }
    
    T y;
    {
      vector<T> v11 = removeIndex(v1, 1);
      vector<T> v22 = removeIndex(v2, 1);
      y = (-1)*det2<T>(v11, v22);
    }
    T z;
    {
      vector<T> v11 = removeIndex(v1, 2);
      vector<T> v22 = removeIndex(v2, 2);
      z = det2<T>(v11, v22);
    }
    return {x, y, z};
  }else{
    printf("Invalid Input.\n");
  }
  return {};
}
 
void printArr2(float* arr, int len, int width){
    vector<vector<float>> vv = arrToVec2(arr, len, width);
    mat m1(vv);
    m1.print();
}
 
template<typename T>
vector<T> subtract(vector<T> v1, vector<T> v2){
  vector<T> v;
  for(int i = 0; i < min(len(v1), len(v2)); i++){
    v.push_back(v1[i] - v2[i]);
  }
  return v;
}
 
template<typename T>
void print(string s, T t){
  fw(s);
  print(t);
}
 
string padRight(string s, string pad, int maxLen){
  int diff = maxLen - len(s);
  int padLen = diff > 0 ? diff : 0;
  string ret = s + repeat(padLen, pad);
  return ret;
}
 
string padLeft(string s, string pad, int maxLen){
  int diff = maxLen - len(s);
  int padLen = diff > 0 ? diff : 0;
  string ret = repeat(padLen, pad) + s;
  return ret;
}
 
//void writeFileRow(string fname, row& r){
//    auto f = [](auto& x) { return toStr(s);};
//    vector<float> vec = {0.1, 1.2};
//    vector<string> s1 = AronLambda::map( f, vec);
//}
 
 
 
 
#endif