//
//  Tensor Regression.cpp
//  Tensor Regression
//
//  Created by Bui  Thang on 7/30/13.
//  Copyright (c) 2013 __MyCompanyName__. All rights reserved.
//

#include <iostream>
#include "Tensor Regression.h"
#include<cstring>
#include<fstream>
#include<string.h>
#include<cmath>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
using namespace std;
int N=0;// The number of siRNA sequence having efficacy score
int N_VH=0,N_Low=0,N_H=0;
int Dim=50;// Maximum number of dimensions
Tensor TensorArray[3000];//Tensor Array to store data set
EncodingMatrix ScoreData[3000];
// The turning paramaters
double Lambda3=0;
Vector alpha,beta,weight;
//-------------------
int n=19;
int K=0;

// Define encoding matrix functions
EncodingMatrix::EncodingMatrix(){ 
    
    for (int i=1;i<=24;i++)
        for (int j=1;j<=4;j++)
            X[i][j]=0;
    label=0;
}
void EncodingMatrix::Get(bool EM[25][5]){
    for (int i=1;i<=n;i++)
        for (int j=1;j<=4;j++)
            EM[i][j]=X[i][j];
}
void EncodingMatrix::Input(string s){
    n=int(s.length());
    for (int i=1;i<=n;i++)
        switch (s[i-1]){
            case 'A': X[i][1]=1; break;
            case 'C': X[i][2]=1; break;
            case 'G': X[i][3]=1; break;
            case 'U': X[i][4]=1; break;
            case 'T': X[i][4]=1; break;
        }
}
void EncodingMatrix::Output(){
    for (int i=1;i<=n;i++){
        for (int j=1;j<=4;j++)
            cout<<X[i][j]<<" ";
        cout<<endl;
    }
}
int EncodingMatrix::Getelement(int p,int q){return X[p][q];}

// Transform scoring siRNA training set to encoding matrices.
void readfile(string filename){
    ifstream obj;
    string str;
    double label;
    int i=1;
    obj.open(filename.c_str());
    while (!obj.eof()) {
        obj>>str;
        obj>>label;
        ScoreData[i++].Input(str);
    }
}
//-----------------------------------
// Define member functions of Transformation matrices class.
Trans_Matrix::Trans_Matrix(){
    
    T=new double*[5];
    for (int i=1; i<=4; i++) {
        T[i]=new double[25];
    }
    for(int i=1;i<=4;i++)
        for (int j=1;j<=24;j++)
            T[i][j]=0;
}
Trans_Matrix::Trans_Matrix(Trans_Matrix &X){
    T=new double*[5];
    for (int i=1; i<=4; i++) {
        T[i]=new double[25];
    }
    for(int i=1;i<=4;i++)
        for (int j=1;j<=24;j++)
            T[i][j]=X.T[i][j];
}
void Trans_Matrix::Set(double a[4][25]){

    for (int i=1; i<=4; i++) {
        for (int j=1; j<=19; j++) {
            T[i][j]=a[i][j];
        }
    }
}
double Trans_Matrix::Getelement(int i,int j){return T[i][j];}
void Trans_Matrix::Initialize(Trans_Matrix Rule){
    double Sum=0;
    double v[5];
    for (int j=1; j<=24; j++){
      //  do{
            Sum=0;
            for(int i=1;i<=4;i++){
                srand (time(NULL));
                T[i][j]= (rand() % 101) * 0.01;
                v[i]=T[i][j];
                Sum+=T[i][j];
            }
            // Standarize T[i][j]
            for(int i=1;i<=4;i++){
                T[i][j]=T[i][j]/Sum;
                v[i]=v[i]/Sum;
            }
     //   }while (!Rule.RuleConstraint(v, j));
    }
}
// Initialize for siRNA design rules.
void Trans_Matrix::Intialize1(){
    double Sum=0;
    for (int j=1; j<=19; j++){
        Sum=0;
        for(int i=1;i<=4;i++){
            srand (time(NULL));
                T[i][j]= (rand() % 101) * 0.001;
                Sum+=T[i][j];

            }
    
        for(int i=1;i<=4;i++)
            T[i][j]=T[i][j]/(Sum);
    }
    Show();
}
Trans_Matrix& Trans_Matrix::operator=(const Trans_Matrix &MT){
    for (int i=1; i<=4; i++) {
        for (int j=1; j<=n; j++) {
            T[i][j]=MT.T[i][j];
        }
    }
    return *this;
}
bool Trans_Matrix::RuleConstraint(double v[],int pos){
    if (pos>19) return true;
    for (int i=1; i<=3; i++) 
        for (int j=i+1; j<=4; j++){
            if(v[i]*v[j]<0) return false;
            if ((T[i][pos]-T[j][pos])!=0)
                if((v[i]<0)&&(v[j]<0)){
                    if (((v[i]-v[j])*(T[i][pos]-T[j][pos]))>0)return false;
                }
                if((v[i]>0)&&(v[j]>0)) if (((v[i]-v[j])*(T[i][pos]-T[j][pos]))<0)return false;
            }
    return true;
}
void Trans_Matrix::Update(double v[],int pos){// Update elements at the column p
    int count=0;
    double sum=0;
    for (int i=1;i<=4;i++){ sum+=v[i];
        if (v[i]<=0) {count++;}
    }

    if ((count==0)||(count==4)){
        for (int i=1;i<=4;i++) T[i][pos]=v[i];
       // cout<<endl;
    }
}
// Show the transformation matrix
void Trans_Matrix::Show(){
    for (int i=1; i<=4; i++) {
        for (int j=1; j<=n; j++) {
            cout<<T[i][j]<<" ";
        }
        cout<<endl;
    }
}
//-------------------------
// Read design rules to transformation matrix;
void RuleRead(string rulefile, Trans_Matrix Rules[]){
    ifstream obj;
    double tmp;
    string str;
    obj.open(rulefile.c_str());
    while (!obj.eof()) {
        ++K;
        for (int i=1; i<=4; i++) 
            for (int j=1; j<=19; j++) {
                obj>>tmp;
                Rules[K].T[i][j]=tmp;
            }
        obj>>str;
    }
}

// Count frequencies of nucleotides at different positions in each class
void FreCounting(string filename, long Fre[5][25],EncodingMatrix Class[], int &num){
    ifstream obj;
    string str;
    int label;
    obj.open(filename.c_str());
   while (!obj.eof()) {
        obj>>str;
        obj>>label;
        num++;
        Class[num].Input(str);
        for (int i=0;i<str.length();i++)
            switch (str[i]) {
                case 'A':
                     ++Fre[1][i+1];
                    break;
                case 'C':
                    ++Fre[2][i+1];
                    break;
                case 'G':
                     ++Fre[3][i+1];
                    break;
                case 'U':
                    ++Fre[4][i+1];
                    break;
                case 'T':
                    ++Fre[4][i+1];
                    break;
                default:
                    break;
            }
    }
}
void VectorDisp(double Al[],int n1){
    for(int i=1;i<=n1;i++) cout<< Al[i]<<"  ";
    cout<<endl;
}
// Display to the screen frequencies of nucleotides at different sequences in each class.
void MatrixDisp(long A[25][25], int n,int m){
    for (int i=1;i<=n;i++){
        for (int j=1; j<=m; j++) {
            cout<<A[i][j]<<" ";
        }
        cout<<endl;
    }
}
double Trans_Matrix::FrobeniusNorm(){
    double Total=0;
    for (int i=1; i<=4; i++) {
        for (int j=1; j<=n; j++) {
            Total+=T[i][j]*T[i][j];
        }
    }
    return sqrt(Total);
}
Trans_Matrix Trans_Matrix::operator-(const Trans_Matrix &MT){
    Trans_Matrix Result;
    for (int i=1; i<=4; i++) {
        for (int j=1; j<=n; j++) {
            Result.T[i][j]=T[i][j]-MT.T[i][j];
        }
    }
    return  Result;
}
// Compute S(p,q); p is the position, q corresponds to nucleotide
//-------------------
// Constructor of Tensor class
Tensor::Tensor(){
    score=0;
    p=0;q=0;
    Arr=NULL;
    Pn=0;
    Prop=NULL;
  
}
Tensor::Tensor(const Tensor &x){
    p=x.p;
    q=x.q;
    Pn=x.Pn;
    score=x.score;
   
    Arr=new double*[p+1];
    for (int i=1; i<=p; i++) {
        Arr[i]=new double[q+1];
    }
    Prop=new double[24];
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Arr[i][j]=x.Arr[i][j];
        }
    }
    for (int i=1;i<=Pn;i++) Prop[i]=x.Prop[i];

}
double Tensor::ProbeniusNorm(){
    double Sum=0;
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Sum+=Arr[i][j]*Arr[i][j];
        }
    }
    return sqrt(Sum);
}
// Build a tensor for the string s with efficacy (value)using K transformation matrices.
void   Tensor::Set(double value, string str, Trans_Matrix TM[]){
    score=double(value);
    double **tmp;
    double *tmp1;
    tmp=Arr;
    int l=p;
    tmp1=Prop;
    
    p=K;
     Pn=0;
    q=int(str.length());
    Arr=new double*[p+1];
    for (int i=1; i<=p; i++) {
        Arr[i]=new double[q+1];
    }
        Prop=new double[Dim];
    for (int i=1; i<=l; i++) {
        delete [] tmp[i];
    }
    delete [] tmp;
    delete [] tmp1;
    for (int i=0; i<str.length(); i++) {
        for(int k=1;k<=K;k++){
            Arr[k][i+1]=0;
        switch (str[i]) {
            case 'A':
                Arr[k][i+1]=TM[k].Getelement(1, i+1);
                break;
            case 'C':
                Arr[k][i+1]=TM[k].Getelement(2, i+1);
                break;
            case 'G':
                Arr[k][i+1]=TM[k].Getelement(3, i+1);
                break;
            case 'U':
                Arr[k][i+1]=TM[k].Getelement(4, i+1);
                break;
            case 'T': Arr[k][i+1]=TM[k].Getelement(4, i+1);
                break;
            default:
                break;
        }
        }
    }
    int GC_Content=0;
    
    for (int i=0; i<str.length(); i++)
        if ((str[i]=='G')||(str[i]=='C')) GC_Content++;
        
    // count GC content fraction
   // Prop[++Pn]=double(GC_Content)/(str.length());
    
    // DeltaT: A/U diffirence between two ends
    int deltaT=0;
    for (int i=0;i<3;i++) {
        if ((str[i]=='A')||(str[i]=='U')) deltaT--;
        if ((str[str.length()-i-1]=='A')||(str[str.length()-i-1]=='U')) deltaT++;
    }
    Prop[++Pn]=double(abs(deltaT))/str.length();

    // check GC tretch constrain on siRNA sequence
    int GC_tretch=0;
    int i=0;
    int j;
    while (i<(str.length())){
        while ((str[i]!='G')&&(str[i]!='C')&&(i<str.length()))i++;
        if (i<str.length()) j=i;
        else break;
        while (((str[j]=='G')||(str[j]=='C'))&&(j<str.length())) j++;
        if (GC_tretch<(j-i)) GC_tretch=j-i;
        i=j;
    }
    //Prop[++Pn]=double(GC_tretch)/str.length();

    // check constrains of A/U contain at 5 end of antisence
    int AU_content1=0;// from position 15 to position 19
    int AU_content2=0;// from position 13 to position 19
    for (int i=14;i<str.length();i++)
        if ((str[i]=='A')||(str[i]=='U')) AU_content1++;
    AU_content2=AU_content1;
    if ((str[12]=='A')||(str[12]=='U')) AU_content2++;
    if ((str[13]=='A')||(str[13]=='U')) AU_content2++;
    //Prop[++Pn]=double(AU_content1)/str.length();
    Prop[++Pn]=double(AU_content2)/str.length();

        //% defaults parameters
        double temp = 25;          //  % temperature in Celsius
        double salt = 0.05;         // % salt concentration in moles per liter (M)
        double primerConc= 50e-6;   // % concentration of primers in mole per liter (M)
       // double weight[5] = {313.21, 289.18, 329.21, 304.2}; //% molecular weight A,C,G,T respectively
        //% compute sequence properties
        int numSeq[20]={0};
        int baseNum[5]={0};

        for (int i=0;i<str.length(); i++)
            switch(str[i]){
                case 'A': numSeq[i+1]=1; baseNum[1]++; break;
                case 'C': numSeq[i+1]=2; baseNum[2]++;break;
                case 'G': numSeq[i+1]=3; baseNum[3]++;break;
                case 'U': numSeq[i+1]=4; baseNum[4]++;break;
                case 'T':numSeq[i+1]=4; baseNum[4]++;break;
                default: break;
            }
       // double w=0;
       // for (int i=1;i<=4;i++) w+= baseNum[i]*weight[i-1];
        
       // w-=61.96;
        double b = 4;
        double basic = 64.9 + (41 * ((baseNum[3] + baseNum[2] - 16.4) / 19)); //TM BASIC [1],[9]
        double saltadj = 100.5 + (41 * (baseNum[3] + baseNum[2])/19.0) - (820.0/19) + (16.6 * log10(salt));// %TM SALT ADJUSTED [9]
        //nearest neighbor parameters from: Panjkovich and Melo, Bioinformatics  Vol 21 no 6 pp 711-722 2004 [1]
        //rows corresponds to A,C,G,T respectively; columns correspond to A,C,G,T respectively
        double Bres86_H[5][5] = {{-9.1,-6.5,-7.8,-8.6},{-5.8,-11,-11.9,-7.8},{-5.6,-11.1,-11,-6.5},{-6,-5.6,-5.8,-9.1}};
        double Bres86_S[5][5] = {{-24,-17.3,-20.8,-23.9},{-12.9,-26.6,-27.8,-20.8},{-13.5,-26.7,-26.6,-17.3},{-16.9,-13.5,-12.9,-24}};
        double Sant96_H[5][5] = {{-8.4,-8.6,-6.1,-6.5},{-7.4,-6.7,-10.1,-6.1},{-7.7,-11.1,-6.7,-8.6},{-6.3,-7.7,-7.4,-8.4}};
        double Sant96_S[5][5] = {{-23.6,-23,-16.1,-18.8},{-19.3,-15.6,-25.5,-16.1},{-20.3,-28.4,-15.6,-23},{-18.5,-20.3,-19.3,-23.6}};
        double Sant98_H[5][5] = {{-7.9,-8.4,-7.8,-7.2},{-8.5,-8,-10.6,-7.8},{-8.2,-9.8,-8,-8.4},{-7.2,-8.2,-8.5,-7.9}};
        double Sant98_S[5][5] = {{-22.2,-22.4,-21,-20.4},{-22.7,-19.9,-27.2,-21},{-22.2,-24.4,-19.9,-22.4},{-21.3,-22.2,-22.7,-22.2}};
        double Sugi96_H[5][5] = {{-8,-9.4,-6.6,-5.6},{-8.2,-10.9,-11.8,-6.6},{-8.8,-10.5,-10.9,-9.4},{-6.6,-8.8,-8.2,-8}};
        double Sugi96_S[5][5] = {{-21.9,-25.5,-16.4,-15.2},{-21,-28.4,-29,-16.4},{-23.5,-26.4,-28.4,-25.5},{-18.4,-23.5,-21,-21.9}};
        
        int ind[19];
        for (int i=1;i<19;i++){
            ind[i]=(numSeq[i+1]-1)*4+numSeq[i];
        }
        // NN is 4x2 matrix. Columns are DeltaH and DeltaS. Rows correspond to
        // methods by Bres86, SantaLucia96, SantaLucia98 and Sugimoto96.
        double NN[5][4]={0};
        int col=0,row=0;
        for (int i=1;i<19;i++){
            col=ind[i]/4;
            row=ind[i]%4;
            if (row!=0)col++;
            if (row==0) row=4;
            col--;row--;
            NN[1][1]+=Bres86_H[row][col];
            NN[1][2]+=Bres86_S[row][col];
            NN[2][1]+=Sant96_H[row][col];
            NN[2][2]+=Sant96_S[row][col];
            NN[3][1]+=Sant98_H[row][col];
            NN[3][2]+=Sant98_S[row][col];
            NN[4][1]+=Sugi96_H[row][col];
            NN[4][2]+=Sugi96_S[row][col];
        }
        
        //% Corrections: all AT pairs, any GC pairs, symmetry, initiation
        //% only AT pairs or any GC pairs?
        bool check=true;
        for(int i=1;i<=19;i++)
            if ((numSeq[i]==2)||(numSeq[i]==3)) {
                check=false;
                break;
            }
        if (check){
            NN[1][2]+=-20.13;
            NN[2][2]+=-9;
            NN[4][1]+=0.6;
            NN[4][2]+=-9;
        }
        else{
            NN[1][2]+=-16.77;
            NN[2][2]+=-5.9;
            NN[4][1]+=0.6;
            NN[4][2]+=-9;
        }
        
        //% initiation with terminal  5'
        
        if((numSeq[1] == 1) || (numSeq[1] == 4)){
            NN[3][1]+=2.3;
            NN[3][2]+=4.1;
        }
        else{
            NN[3][1]+=0.1;
            NN[3][2]+=-2.8;
            
        }
        
        if ((numSeq[19]==2)||(numSeq[19]==3)){
            NN[3][1]+=0.1;
            NN[3][2]+=-2.8;
        }
        else{
            NN[3][1]+=2.3;
            NN[3][2]+=4.1;
        }
        
        double tm[7]={0};
        tm[1]=basic;
        tm[2]=saltadj;
        for (int i=3;i<=6;i++)
            tm[i] = (NN[i-2][1] * 1000.0/(NN[i-2][2] + (1.9872 * log(primerConc/b)))) + (16.6 * log10(salt)) - 273.15; //%TM NEAREST NEIGHBOR
        
        for (int i=1;i<=4;i++)
            NN[i][3] = NN[i][1] - ((temp+273.15)* (NN[i][2]/1000));// % DeltaG
        // show result;
       // Prop[++Pn]=w/10000;
        
        // Melting temperature properties
        for (int i=1;i<=4;i++)Prop[++Pn]=tm[i]/1000;//i=1..6 (i:1..4 good)
        // Thermodynamic properties 
       // for (int i=1;i<=4;i++)
        for (int j=1;j<=2;j++) Prop[++Pn]=NN[1][j]/1000;//j1=-->3
   
 
}
// Get a tensor
void Tensor::Get(double Matrix[][25],int &n,int &m){
    n=p;
    m=q;
    for (int i=1; i<=n; i++) {
        for (int j=1; j<=m; j++) {
            Matrix[i][j]=Arr[i][j];
        }
    }
}
// Consider a array as tensor;
void Tensor::Input(double **Matrix,int n, int m){
    double **tmp;
    tmp=Arr;
    int l=p;
        p=n;
        q=m;
        Arr=new double*[p+1];
        for (int i=1; i<=p; i++) {
            Arr[i]=new double[q+1];
        }
    for (int i=1; i<=l; i++) {
        delete [] tmp[i];
    }
   
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
           Arr[i][j]= Matrix[i][j];
        }
    }   
}
// Display tensor
ostream &operator<<(ostream &out, const Tensor x ){
    out<<x.score<<endl;
    for (int i=1; i<=x.p; i++) {
        for (int j=1; j<=x.q; j++) {
            out<< x.Arr[i][j]<<" ";
        }
        out<<endl;
    }
    return  out;
}
// Transform siRNA sequences in data set in to Tensor and store in TensorArray
void TensorData(string rulefile, Trans_Matrix TM[],Tensor Dataset[],int &L, char sign){
    ifstream obj;
    string str;
    double tmp;
    L=0;
    obj.open(rulefile.c_str());
    tmp=0;
    while (!obj.eof()) {
        ++L;
        obj>>str;
        if (sign!='0')obj>>tmp;
        fflush(stdin);
       // cout<<str<<endl;
        if (str!="") Dataset[L].Set(tmp, str, TM);
    }
    L--;
    obj.close();
}
// Display tensor data set
void Display(string rulefile, Trans_Matrix TM[]){
    TensorData(rulefile, TM,TensorArray,N,'0');
    for (int i=1; i<=N; i++) {
        cout<<TensorArray[i];
        system( "read -n 1 -s -p \"Press any key to continue...\"" );
        cout<<endl;
    }
}
// Get value of the tensor at the cell (i,j)
double Tensor::value(int i,int j)const{return Arr[i][j];}
// Get x axis
int Tensor::Getx()const{return p;};
// Get y axis
int Tensor::Gety()const{return q;};
// Overloading operator=
Tensor &Tensor::operator=(const Tensor &y){
   
     if (Arr==NULL) {
        p=y.p;
        q=y.q;
         Pn=y.Pn;
        Arr=new double*[p+1];
        for (int i=1; i<=p; i++) {
        Arr[i]=new double[q+1];
        }
         Prop= new double[Pn+1];
         
    }
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Arr[i][j]=y.Arr[i][j];
        }
    }
    for (int i=1; i<=Pn; i++) {
        Prop[i]=y.Prop[i];
    }
    return *this;
}
Tensor Tensor::operator=(const double a){
    for (int i=1; i<=p; i++)
        for (int j=1;j<=q;j++)
        Arr[i][j]=a;
    
    return *this;
    
};
// Get efficacy score
double Tensor::GetScore()const{return score;}
// Overloading operator+: Tensor + constant
Tensor Tensor::operator+(double scale)const{
    Tensor Add;
    Add.p=p;
    Add.q=q;
    Add.Arr=new double*[p+1];
    for (int i=1; i<=p; i++) {
        Add.Arr[i]=new double[q+1];
    }
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Add.Arr[i][j]=Arr[i][j]+scale;
        }
    }
    return Add;
}
void Tensor::operator+=(const double scale){
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Arr[i][j]=Arr[i][j]+scale;
        }
    }

}
void Tensor::operator+=(const Tensor &y){
    
    if (Arr==NULL){
        p=y.p;
        q=y.q;
        Arr=new double*[p+1];
        for (int i=1; i<=p; i++) {
            Arr[i]=new double[q+1];
        }
        for (int i=1; i<=p; i++)
            for (int j=1; j<=q; j++)
                Arr[i][j]=y.Arr[i][j];
            
    }
     else {
            for (int i=1; i<=p; i++) {
            for (int j=1; j<=q; j++) {
                Arr[i][j]=Arr[i][j]+y.Arr[i][j];
            }
        }
     }
}
Tensor operator+(const Tensor &x, const Tensor &y){
    Tensor result;
    result.p=y.p;
    result.q=y.q;
    result.Arr=new double*[y.p+1];
    for (int i=1; i<=y.p; i++) {
        result.Arr[i]=new double[y.q+1];
    }

    for (int i=1; i<=y.p; i++)
        for (int j=1; j<=y.q; j++)
            result.Arr[i][j]=x.Arr[i][j]+y.Arr[i][j];
    

    return result;
};
Tensor Tensor::operator*(double scale)const{
    Tensor Product;
    Product.p=p;
    Product.q=q;
    Product.Arr=new double*[p+1];
    for (int i=1; i<=p; i++) {
        Product.Arr[i]=new double[q+1];
    }
    for (int i=1; i<=p; i++) {
        for (int j=1; j<=q; j++) {
            Product.Arr[i][j]=Arr[i][j]*scale;
        }
    }
    return Product;
}
Tensor Tensor::Inverse()const{
    double ratio,a;
    double **matrix=NULL;
    bool ok=true;
    matrix=new double*[2*p+1];
    for (int i=1;i<=2*p;i++) matrix[i]=new double[2*p+1];
    
    for(int i=1;i<=p;i++)
        for (int j=1;j<=p;j++) {
            matrix[i][j]=this->value(i,j);
            if (matrix[i][j]!=0) ok=false;
    }
    
    for(int i = 1; i <= p; i++){
        for(int j = p+1; j <= 2*p; j++){
            if(i==(j-p))
                matrix[i][j] = 1.0;
            else
                matrix[i][j] = 0.0;
        }
    }
    for(int i = 1; i <=p; i++){
        for(int j = 1; j <=p; j++){
            if(i!=j){
                a=matrix[i][i];
                if (a==0) {
                    a=0.0001;
                }
               
                ratio = matrix[j][i]/a;//(matrix[i][i]);
                for(int k = 1; k <=2*p; k++){
                    matrix[j][k] -= ratio * matrix[i][k];
                }
            }
        }
    }
    for(int i = 1; i <= p; i++){
        a = matrix[i][i];
        if (a==0) {
            a=0.0001;
        }
        for(int j = 1; j <=2*p; j++){
            matrix[i][j] /= (a);
        }
    }
    double **InvMatrix=NULL;
    InvMatrix=new double*[p+1];
    for (int i=1;i<=p;i++) InvMatrix[i]=new double[p+1];
    if(ok){
        for(int i = 1; i<=p; i++)
            for(int j = p+1; j <=2*p; j++)
                InvMatrix[i][j-p]=0;
    }
    else{
    for(int i = 1; i<=p; i++)
        for(int j = p+1; j <=2*p; j++)
            InvMatrix[i][j-p]=matrix[i][j];
    }
    
    Tensor x;
    x.Input(InvMatrix, p, p);
 for (int i=1;i<=2*p;i++) delete[] matrix[i];
   for (int i=1;i<=p;i++) delete[] InvMatrix[i];
    return x;
}
Vector Tensor::GetP(){
    Vector ve;
    ve.Set(Prop, Pn);
    return ve;
}
// Constructor of Vector class
Vector::Vector(){
    len=0;
    v=new double[Dim];
    for (int i=1; i<Dim; i++) v[i]=0;
}
Vector::Vector(const Vector &x){
    len=x.len;
    v=new double[len+1];
    for (int i=1; i<=len; i++) {
        v[i]=x.v[i];
    }
}
Vector::~Vector(){delete []v;}
// Initialize a vector
void Vector::Initialize(int x){
    len=x;
    for (int i=1; i<=len; i++) {
        srand (time(NULL));
        v[i]= (rand() % 101) * 0.1;
    }   

}
void Vector::Set(double a[],int l){
    len=l;
    for (int i=1; i<=l; i++) {
        v[i]=a[i];
    }
}
// Product of 2 vectors
Tensor Vector::operator*(const Vector y){
    double **Result=NULL;
    Tensor x;
    Result=new double*[len+1];
    for (int i=1; i<=len; i++) {
        Result[i]=new double[y.len+1];
    }
    for (int i=1; i<=len; i++) {
        for (int j=1; j<=y.len; j++) {
            Result[i][j]=v[i]*y.v[j];
        }
    }
    x.Input(Result, len, y.len);
    for (int i=1;i<=len;i++) delete[] Result[i];
    return x;
}
// Product of a tensor and vector
Vector operator*(const Tensor &y,const Vector &x){
    Vector z;
    z.len=y.Getx();
    for (int i=1; i<=z.len; i++){
        z.v[i]=0;
        for (int j=1; j<=x.len; j++) 
           z.v[i]+=x.v[j]*y.value(i,j);
    }
    return z;
}
// Product of a vector and a tensor
Vector operator*(const Vector &x, const Tensor &y){
    Vector z;
    z.len=y.Gety();
    for (int j=1; j<=z.len; j++){
        z.v[j]=0;
        for (int i=1; i<=y.Getx(); i++) 
            z.v[j]+=x.v[i]*y.value(i,j);
    }
    return z;
}
// Scale vector with constant y
Vector Vector::operator*(double y)const{
    Vector z;
    z.len=len;
    for (int j=1; j<=len; j++) 
            z.v[j]=v[j]*y;
    return z;
}
void Vector::operator*=(const double y){
    for (int j=1; j<=len; j++)
        v[j]=v[j]*y;
};
double Vector::InnerProduct(const Vector &y){
    double Sum=0;
    for (int i=1; i<=len; i++) {
        Sum+=v[i]*y.v[i];
    }
    return Sum;
}
// Subtraction of 2 vector
Vector Vector::operator-(const Vector &y)const{
    Vector z;
    z.len=len;
    for (int i=1; i<=len; i++) {
        z.v[i]=v[i]-y.v[i];
    }
    return z;
}
// Sum of 2 vectors
Vector Vector::operator+(const Vector &y)const{
    Vector z;
    z.len=len;
        for (int i=1; i<=len; i++) {
            z.v[i]=v[i]+y.v[i];
        }
    return z;
}
void Vector::operator+=(const Vector &y){
 
    len=y.len;
    for (int i=1; i<=len; i++) {
        v[i]=v[i]+y.v[i];
    }
    
}
void Vector::operator=(const Vector &y){
    len=y.len;
    for (int i=1; i<=y.len; i++) {
        v[i]=y.v[i];
    }
    
}
void Vector::operator=(const double a){
    for(int i=1;i<=len;i++) v[i]=a;
   
    
}
void Vector::Normalize(){
    for (int i=1; i<=len; i++) {
        if(abs(v[i])<pow(10.0, -5)) v[i]=0;
    }
}
void Vector::standlize(){
    double Sum=0;
    for (int i=1; i<=len; i++) {
        Sum+=v[i];
    }
    for (int i=1; i<=len; i++) v[i]=v[i]/Sum;
    
    
}
ostream &operator<<(ostream &out,const Vector y){
    for (int i=1; i<=y.len; i++) {
        out<<y.v[i]<<"  ";
    }
    return  out;
}
double Vector::Norm2(){
    double S=0;
    for (int i=1; i<=len; i++) {
        S+=v[i]*v[i];
    }
    return sqrt(S);
}

// Partition dataset into training set and testing set: Ten is the number of testing set, Trn is the number of training set
void Gettest(int indices[], int Test[], int &num, int k, int Train[],int &Trn){
    num=Trn=0;
    for (int i=1; i<=N; i++) 
        if (indices[i]==k) {
            num++;
            Test[num]=i;
            //   cout<<i<<" ";
        }
        else{Train[++Trn]=i;
        }
}
// Tensor Regression predictor;
void TensorPredictor(int Test[],int Ten,Vector &alpha, Vector &beta, double Target[], double Pred[]){
    
    for(int i=1;i<=Ten;i++){
        Target[i]=TensorArray[Test[i]].GetScore();
        Pred[i]=(alpha*TensorArray[Test[i]]).InnerProduct(beta)+Lambda3*weight.InnerProduct(TensorArray[Test[i]].GetP());
    }
    /*for (int i=1; i<=Ten; i++) {
     cout<<Target[i]<<"   "<<Pred[i]<<endl;
     }*/
}


void Predictor(Tensor Testset[],int Ntest,Vector &alpha, Vector &beta, double Target[], double Pred[]){
    
    for(int i=1;i<=Ntest;i++){
        Target[i]=Testset[i].GetScore();
        Pred[i]=(alpha*Testset[i]).InnerProduct(beta)+Lambda3*weight.InnerProduct(Testset[i].GetP());
    }
}
// Correlation Coefficient of 2 vectors;
double Corr(double x[],double y[], int size){
    double Mx=0,My=0;
    for (int i=1; i<=size; i++) {
        Mx+=x[i];
        My+=y[i];
    }
    Mx=Mx/size;
    My=My/size;
    double numerator=0, denominator1=0,denominator2=0;
    for (int i=1; i<=size; i++) {
        numerator+=(x[i]-Mx)*(y[i]-My);
        denominator1+=(x[i]-Mx)*(x[i]-Mx);
        denominator2+=(y[i]-My)*(y[i]-My);
    }
    return numerator/(sqrt(denominator1)*sqrt(denominator2));
}

void Getdata(string input,string output, string label, Trans_Matrix &R){
    ifstream obj;
    string s;
    ofstream out,lab;
    int tmp=0;
    obj.open(input.c_str());
    out.open(output.c_str(),ios::app);
    lab.open(label.c_str(),ios::app);
    while (!obj.eof()) {
        obj>>tmp;
        obj>>s;
       // cout<<s<<endl;
        lab<<tmp<<endl;
        for (int i=0; i<s.length(); i++)
        switch (s[i]) {
            case 'A':
                out<<R.Getelement(1, i+1)<<" ";
                break;
            case 'C':
                out<<R.Getelement(2, i+1)<<" ";
                break;
            case 'G':
                out<<R.Getelement(3, i+1)<<" ";
                break;
            case 'U':
                out<<R.Getelement(4, i+1)<<" ";
                break;
            case 'T':
                out<<R.Getelement(4, i+1)<<" ";
                break;
            default:
                break;
        }
        out<<endl;
    }
    obj.close();
    out.close();
    lab.close();
}
void writetofile(Trans_Matrix T[]){

    ofstream obj;
    obj.open("Output.strings",ios::app);
   // obj<<"Transformation Matrices:"<<endl;
    for (int k=1; k<=K; k++) {
        for (int i=1; i<=4; i++) {
            for (int j=1; j<=19; j++) {
                obj<<T[k].Getelement(i,j)<<" ";
            }
            obj<<endl;
        }
//        obj<<endl;
    }
  //  obj<<"Alpha:"<<endl;
    for (int k=1; k<=K;k++) {
        obj<<alpha[k]<<" ";
    }
    obj<<endl;
   // obj<<endl<<"beta:"<<endl;
    for (int k=1; k<=n;k++) {
        obj<<beta[k]<<" ";
    }
    obj<<endl;
   // obj<<endl<<"Weight:"<<endl;
    for (int k=1; k<=weight.length();k++) {
        obj<<weight[k]<<" ";
    }
    obj<<endl<<Lambda3<<endl;
    obj.close();
    
}
void Test(string testfile, char sign){
    Trans_Matrix T[50];
    K=7;
    Tensor DTest[3000];//,HUTest[500],HaBoTest[500],ReyTest[500],VicTest[500];//,Uitei[200],Sloan[200];
    int Ntest=0;//HUn=0,HAn=0,VIn=0,REn=0;//,Uin=0,Sn=0;
    double Target[3000],Pred[3000];
    ifstream obj;
    string str;
    double a[5][25];
    double tmp[25];
    int size;
    double Correlation=0;
  //  double t1=0,t2=0,t3=0;//,t4=0,t5=0;
    obj.open("ParaModel.strings");
    alpha.Initialize(K);
    int l=0;
    while (!obj.eof()) {
        l++;
       // obj>>str;
               for (int k=1; k<=K; k++) {
            for (int i=1; i<=4; i++) {
                for (int j=1; j<=19; j++) {
                    obj>>a[i][j];
                }
            }
            T[k].Set(a);
         //   T[k].Show();
          //  cout<<endl;
        }
      //  system( "read -n 1 -s -p \"Press any key to continue...\"" );
        
        for (int k=1; k<=K; k++) {
            obj>>tmp[k];
        }
        alpha.Set(tmp, K);
       // cout<<alpha<<endl;
        for (int j=1; j<=n; j++) {
            obj>>tmp[j];
            
        }
        fflush(stdin);
        beta.Set(tmp, n);
      //  cout<<"Beta:"<<beta<<endl;
        for (int i=1; i<=8; i++) {
            obj>>tmp[i];
        }
        size=8;
        weight.Set(tmp,size);
      //  cout<<"Weight="<<weight<<endl;
        obj>>Lambda3;
       // cout<<Lambda3<<endl;
        //if (l==k	) {// 66
            
        
        TensorData(testfile,T,DTest,Ntest,sign);
        Predictor(DTest,Ntest, alpha, beta, Target, Pred);
        //cout<<Ntest<<endl;
        if (sign=='0')
        for (int i=1; i<=Ntest; i++) {
            cout<<Pred[i]<<endl;
        }
            
        else {
           // cout<<testfile<<endl;
            Correlation=Corr(Target, Pred, Ntest);
            
            cout<<Correlation<<"  \t"<<endl;
           // cout<<"Solution "<<l<<endl;
            for (int i=1; i<=Ntest; i++) {
                cout<<Target[i]<<"\t"<<Pred[i]<<endl;
            }
           // if (l==85)system( "read -n 1 -s -p \"Press any key to continue...\"" );
        }
      // }
        //system( "read -n 1 -s -p \"Press any key to continue...\"" );
       
       // system( "read -n 1 -s -p \"Press any key to continue...\"" );

        //if (t>1) check=false;
        
    }
   // cout<<l<<endl;
    obj.close();

}