void DNN::backPropSet(const dvec& input,const dvec& output)
{

    unsigned int i;

    feedForward(input);

    unsigned int L = activations.size() - 1;

    /*
     * Start with the final layer
     */
    std::vector<Matrix> d(L+1);

    /*
     * Copy contents
     */
    Matrix out(output.size(),1);

    for(i=0;i<output.size();++i) out(i,0) = output.at(i);

    /*
     * Final layer error
     */
    Matrix DC = Matrix::apply(quad_dC_dA,activations.at(L),out);
    d.at(L)   = Matrix::had(DC,activations.at(L));

    /*
     * Backpropagate
     */
    for(i=L;i>0;--i)
    {

        Matrix wd = weights.at(i-1).T() * d.at(i);
        d.at(i-1)    = Matrix::had( wd, activations.at(i-1) );

    }

    /*
     * Calculate the gradient cost for this set
     */
    for(i=L;i>0;--i)
    {

        bGradient.at(i-1) = bGradient.at(i-1) + d.at(i);

        Matrix wg = d.at(i) * activations.at(i-1).T();
        wGradient.at(i-1) = wGradient.at(i-1) + wg;

    }

}

int gauss(dvec& x, dvec& w) {

    int n = x.size();
    double dist = 1;
    double tol = 1e-15;
    int iter = 0;

    // Use Chebyshev-Gauss nodes and initial guess
    initguess(x);
//    chebnodes(x);

    dvec x0(x);
    dvec L(n,0.0);
    dvec Lp(n,0.0);
    dvec a(n,0.0);
    dvec b(n,0.0);

    rec_legendre(a,b);

    // Iteratively correct nodes with Newton's method
    while(dist>tol) {     
        newton_step(a,b,x,L,Lp);
        dist = dist_max(x,x0); 
        ++iter;
        x0 = x;
    } 

    // Compute Weights
    for(int i=0;i<n;++i){
        w[i] = 2.0/((1-x[i]*x[i])*(Lp[i]*Lp[i]));
    }

    return iter;
}

// Update grid points as well as the nth Legendre polynomial and its derivative  
void newton_step(const dvec& a, const dvec& b, dvec& x, dvec& L, dvec& Lp) {
    int n = x.size();
    legendre(a,b,x,L,Lp);
    for(int i=0;i<n;++i) {
        x[i] -= L[i]/Lp[i];
    } 
}

// Compute Legendre polynomial recursion coefficients
void rec_legendre(dvec& a, dvec& b){
    int n = a.size();
    for(int i=0;i<n;++i) {
        a[i] = (2.0*i+1.0)/(i+1.0);
        b[i] = double(i)/(i+1.0);
    }
}

// Compute Chebyshev Gauss Nodes
double chebnodes(dvec& x) {
    int n = x.size();

    for(int i=0;i<n;++i) {
        x[i] = cos(pi*double(2*i+1)/double(2*n));    
    } 
} 

	void backProp(int target, Delta & d) {
		dvec<OUTPUT_SIZE> dErr;
		//derive loss function at softmax
		int ans = 0;
		for (int i = 0; i < OUTPUT_SIZE; ++i) {
			dErr.d[i] = output.d[i];
			if (output.d[i] > output.d[ans]) {
				ans = i;
			}
			if (target == i) dErr.d[i] -= 1;
		}
		//derive at hidden
		dvec<HIDDEN_SIZE> dHidden;
		MVLeftMultiply(model->Wfc, dErr, dHidden);
		//derive at output gate
		dvec<HIDDEN_SIZE> dGo(dHidden);
		dGo *= cellTanh; dGo *= Go.getDerivSigmoid();
		//derive at cell
		dvec<HIDDEN_SIZE> dCell(dHidden);
		dCell *= Go; dCell *= cellTanh.getDerivTanh();
		//derive at input/forget/hidden gate
		dvec<HIDDEN_SIZE> dGi(dCell); dGi *= Gu; dGi *= Gi.getDerivSigmoid();
		dvec<HIDDEN_SIZE> dGu(dCell); dGu *= Gi; dGu *= Gu.getDerivTanh();
		dvec<HIDDEN_SIZE> dGf(dCell); dGf *= lastState->cell; dGf *= Gf.getDerivSigmoid();

		//derive at input embedding
		d.input.initToZero();
		dvec<EMBEDDING_SIZE> tmp;
		MVLeftMultiply(model->Wi, dGi, tmp); d.input += tmp;
		MVLeftMultiply(model->Wf, dGf, tmp); d.input += tmp;
		MVLeftMultiply(model->Wo, dGo, tmp); d.input += tmp;
		MVLeftMultiply(model->Wu, dGu, tmp); d.input += tmp;
		//derive for bias
		d.Bf = dGf; d.Bi = dGi;
		d.Bo = dGo; d.Bu = dGu;

		//derive for weights
		VVOuterProduct(dGo, *input, d.Wo);
		VVOuterProduct(dGi, *input, d.Wi);
		VVOuterProduct(dGf, *input, d.Wf);
		VVOuterProduct(dGu, *input, d.Wu);
		VVOuterProduct(dGo, lastState->hidden, d.Uo);
		VVOuterProduct(dGi, lastState->hidden, d.Ui);
		VVOuterProduct(dGf, lastState->hidden, d.Uf);
		VVOuterProduct(dGu, lastState->hidden, d.Uu);
	}

void DNN::feedForward(const dvec & inputs)
{

    unsigned int layers = activations.size();
    unsigned int i;

    for(i=0;i<inputs.size();++i)
    {
        activations.at(0)(i,0) = sigmoid(inputs.at(i));
    }

    for(i=1;i<layers;++i)
    {
        activations.at(i) = (weights.at(i-1)*activations.at(i-1) + bias.at(i-1)).apply(sigmoid);
    }

}

double initguess(dvec& x) {
    int n = x.size();

    for(int i=0;i<n;++i) {
        x[i] = -cos(pi*double(i+.75)/double(n+0.5));    
    } 

    
}

double simple_energy( const dvec &q , const dvec &p )
{
    using boost::math::pow;
    const size_t N=q.size();
    double e(0.0);
    for( size_t n=0 ; n<N ; ++n )
    {
        e += 0.5*pow<2>( p[0] ) 
            + pow<4>( q[0] )/4.0;
    }
    return e;
}

	void feedForward(LSTMState *lastState, const string& inputStr) {
		this->lastState = lastState;
		//get input char embedding
		input = &model->lookupTable[inputStr];
		
		//lstm
		dvec<HIDDEN_SIZE> tmp;
		
		MVRightMultiply(model->Wi, *input, tmp);
		MVRightMultiply(model->Ui, lastState->hidden, Gi);
		Gi += tmp; Gi += model->Bi;  Gi.performSigmoid();
		
		MVRightMultiply(model->Wf, *input, tmp);
		MVRightMultiply(model->Uf, lastState->hidden, Gf);
		Gf += tmp; Gf += model->Bf; Gf.performSigmoid();

		MVRightMultiply(model->Wo, *input, tmp);
		MVRightMultiply(model->Uo, lastState->hidden, Go);
		Go += tmp; Go += model->Bo; Go.performSigmoid();

		MVRightMultiply(model->Wu, *input, tmp);
		MVRightMultiply(model->Uu, lastState->hidden, Gu);
		Gu += tmp; Gu += model->Bu; Gu.performTanh();

		cell = lastState->cell; cell *= Gf;
		tmp = Gi; tmp *= Gu;
		cell += tmp;

		cellTanh = cell; cellTanh.performTanh();
		hidden = cellTanh; hidden *= Go;

		//output to fc layer
		MVRightMultiply(model->Wfc, hidden, output);
		output.performSoftmax();
	}

// Compute maximum pointwise difference
double dist_max(const dvec& a, const dvec& b){
    int n = a.size();
    double dist;
    double max_dist;

    for(int i=0; i<n; ++i) {
        dist = std::abs(a[i]-b[i]);
        if(dist>max_dist) {
            max_dist = dist; 
        }  
    }
    return max_dist;
}

 dvec operator()( const dvec &x , dvec dxdt )
 {
     const size_t N = x.size();
     //hpx::cout << boost::format("rhs start %d , %d \n") % x.size() % dxdt.size() << hpx::flush;
     dxdt[0] = coupling( x[1]-x[0] );
     for( size_t i=1 ; i<N-1 ; ++i )
     {
         dxdt[i] = coupling( x[i+1]-x[i] ) + coupling( x[i-1]-x[i] );
     }
     dxdt[N-1] = coupling( x[N-2] - x[N-1] );
     //hpx::cout << "rhs end\n" << hpx::flush;
     return dxdt;
 }

dvec DNN::predict(const dvec & inputs)
{

    if(inputs.size() != activations.at(0).rows())
        throw std::runtime_error("Input size must equal the number of input neurons");

    feedForward(inputs);

    unsigned int layers = activations.size();

    return activations[layers-1].col_slice(0,0,activations[layers-1].rows()-1);

}

void ConvNet::learn(dmatrix3 &stimulus, dvec &target)
{
    dvec result(target.size());
    
    result = feedforward(stimulus);
    real er = 0;
    for(int x=0;x<target.size();x++) {
        L4.Errors[x] = sigmoid_p(target[x]) * (target[x] - result[x]);
        er += L4.Errors[x];
    }
    
    std::cout << "Output error: " << er << std::endl;
    
    L3.Errors = L4.backpropagation();
    dvec l2er = L3.backpropagation();
    L2.Errors = fold3(l2er, L2.OutShape);
    L1.Errors = L2.backpropagation();
    
    L4.weight_update(L3.Activations);
    L3.weight_update(flatten(L2.Activations));
    //L2 is a PoolLayer, those doesn't have any weights.
    L1.weight_update(Inputs); // Warning: Input is a pointer!!!
}

// Evaluate the nth order Legendre polynomial and its derivative
void legendre(const dvec& a, const dvec& b, const dvec& x, dvec& L, dvec& Lp) {
    int n = x.size();
    dvec L0(n,1.0);
    dvec L1(n,0.0);
    
    // Iterate over grid points
    for(int j=0;j<n;++j) {
        L1[j] = x[j];
        // Iterate over polynomials  
        for(int k=1;k<n;++k) {

            L[j] = a[k]*x[j]*L1[j]-b[k]*L0[j];
            L0[j] = L1[j];
            L1[j] = L[j];
        }
        Lp[j] = n*(L0[j]-x[j]*L[j])/(1.0-x[j]*x[j]); 
    } 
} 

double energy( const dvec &q , const dvec &p )
{
    using std::pow;
    using std::abs;
    using boost::math::pow;
    const size_t N=q.size();
    double e(0.0);
    for( size_t n=0 ; n<N-1 ; ++n )
    {
        e += 0.5*pow<2>( p[n] ) 
            + pow( abs(q[n]) , KAPPA )/KAPPA 
            + pow( abs(q[n]-q[n+1]) , LAMBDA )/LAMBDA;
    }
    e += 0.5*pow<2>( p[N-1] ) 
        + pow( abs(q[N-1]) , KAPPA )/KAPPA 
        + pow( abs(q[N-1]-q[0]) , LAMBDA )/LAMBDA;
    return e;
}

void CanteraGas::getDestructionRates(dvec& wDot) const
{
    kinetics->getDestructionRates(wDot.data());
}

void CanteraGas::getReactionRates(dvec& wDot) const
{
    kinetics->getNetProductionRates(wDot.data());
}

void CanteraGas::getEnthalpies(dvec& hk) const
{
    thermo.getPartialMolarEnthalpies(hk.data());
}

void CanteraGas::getSpecificHeatCapacities(dvec& cpSpec) const
{
    thermo.getPartialMolarCp(cpSpec.data());
}