$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2011 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#ifndef QUB_CONSTRAINTS_H
#define QUB_CONSTRAINTS_H

#include "ublas_plus.h"

#include "qubfast.h"

#include "callbk_reportfun.h"

/*

*/


/*
 What are the optimization parameters?  (some of) the off-diagonal elements of K0 and K1.
 Actually, we use log(K0) to keep them positive.
 Breaking with past practice, we include all the unconnected zeros,
 and put all the K0s before all the K1s.  As we'll see, this is a column vector.
*/

extern "C" QUBFAST_API int CountK(int Ns);
extern "C" QUBFAST_API void K01toK(int Ns, double **K0, double **K1, double *K);
extern "C" QUBFAST_API void K012toK(int Ns, double **K0, double **K1, double **K2, double *K);

// K alone isn't enough to reconstruct K0 and K1, since 0 could mean log(1), or else unconnected.
// To disambiguate, you provide a valid K0 to be updated -- its zeros mean unconnected.

extern "C" QUBFAST_API int KtoK01(int Ns, double *K, double **K0, double **K1);
extern "C" QUBFAST_API int KtoK012(int Ns, double *K, double **K0, double **K1, double **K2);
// returns 0 on success, or -(i+1) if exp(K[i]) is too tiny or huge

extern "C" QUBFAST_API int dK_to_dK01(int Ns, double *K, double *dK, double **K0, double **d_K0, double **d_K1);
extern "C" QUBFAST_API int dK_to_dK012(int Ns, double *K, double *dK, double **K0, double **d_K0, double **d_K1, double **d_K2);

// But some of those K are zero and have to stay zero, and the user may have other constraints.
// We represent the constraints by the system of linear equations (11) and derive the transformations (12, 13).

// Ncns_max = upper_bound(Ncns) = Nk [fixed] + 2*Ncns_user
// Ain[Ncns_max][Nk]
// Bin[Ncns_max]
// Ncns = 0 (or how many pre-existing constraints in Ain,Bin)

// Some constraints are more a part of the representation than the model:
//   * where K0 == 0 there's no connection; fix k0 and k1
//   * where V == 0  fix k1

extern "C" QUBFAST_API void AddAutoConstraints(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
							int *Ncns, double **Ain, double *Bin);

// The user might have some constraints too:

extern "C" QUBFAST_API void AddFixRateConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						 int *Ncns, double **Ain, double *Bin,
						 int a, int b);
extern "C" QUBFAST_API void AddFixExpConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V,int **P,
						int *Ncns, double **Ain, double *Bin,
						int a, int b);
extern "C" QUBFAST_API void AddFixPresConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V,int **P,
						 int *Ncns, double **Ain, double *Bin,
						 int a, int b);
extern "C" QUBFAST_API void AddScaleRateConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						   int *Ncns, double **Ain, double *Bin,
						   int a, int b, int c, int d);
extern "C" QUBFAST_API void AddScaleExpConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						  int *Ncns, double **Ain, double *Bin,
						  int a, int b, int c, int d);
extern "C" QUBFAST_API void AddScalePresConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						   int *Ncns, double **Ain, double *Bin,
						   int a, int b, int c, int d);
extern "C" QUBFAST_API int AddLoopBalConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						int *Ncns, double **Ain, double *Bin,
						int Nloop, int *loopstates,
						callbk_reportfun report_cb, void *report_data);
extern "C" QUBFAST_API int AddLoopImbalConstraint(int Ns, double **K0, double **K1, double **K2, int **L, int **V, int **P,
						  int *Ncns, double **Ain, double *Bin,
						  int Nloop, int *loopstates,
						  callbk_reportfun report_cb, void *report_data);

// It is a lot faster if you first eliminate fixed parameters which don't participate in any other constraints:
extern "C" QUBFAST_API void IdentifyFixedParams_p(int Nk, int Ncns, double *K, double *Ain, double *Bin,
			       	/*output: */   	int *NkR, int *kR_indices, int *NcnsR, int *constraint_indices);
extern "C" QUBFAST_API void IdentifyFixedParams(int Nk, int Ncns, double *K, double **Ain, double *Bin,
			       	/*output: */   	int *NkR, int *kR_indices, int *NcnsR, int *constraint_indices);
extern "C" QUBFAST_API void EliminateFixedParams_p(int Nk, int Ncns, double *K, double *Ain, double *Bin,
						   int NkR, int *kR_indices, int NcnsR, int *constraint_indices,
						   double *Kr, double *AinR, double *BinR);
extern "C" QUBFAST_API void EliminateFixedParams(int Nk, int Ncns, double *K, double **Ain, double *Bin,
						 int NkR, int *kR_indices, int NcnsR, int *constraint_indices,
						 double *Kr, double **AinR, double *BinR);
extern "C" QUBFAST_API void Kr_to_K(double *Kr, int NkR, double *K, int *kR_indices);

/*
 SetupLinearConstraints requires that all constraints be linearly independent.
 First call ReduceConstraints to delete any dependent constraints.
*/
extern "C" QUBFAST_API void ReduceConstraints(int Nk, int *Ncns, double **Ain, double *Bin);
extern "C" QUBFAST_API int ReduceConstraints_p(int Nk, int Ncns, double *Ain, double *Bin);


/*
 Now that we have K0, K1, L, V and rate vector K,  
 we do some magic with the SVD to change your linear constraints
    A * K = B
 into a transformation between K and the column vector P of free parameters
    K = A*P + B
    P = Ainv * K
*/

extern "C" QUBFAST_API int SetupLinearConstraints(int Nc, double **Ain, double *Bin,
							    int Nk, double *K, double *P_out, // Nk = Nc + Np
							    double **Acns, double *Bcns, double **Ainv,
								callbk_reportfun report_cb, void *report_data);
// Acns[Nk][Np]; Bcns[Nk]; Ainv[Np][Nk]

extern "C" QUBFAST_API int SetupLinearConstraints_p(int Nc, double *Ain, double *Bin,
							      int Nk, double *K, double *P_out, // Nk = Nc + Np
							      double *Acns, double *Bcns, double *Ainv,
							      callbk_reportfun report_cb, void *report_data);


extern "C" QUBFAST_API void KtoPars(int Nk, int Np, double *K, double **A, double **Ainv, double *B, double *pars);

extern "C" QUBFAST_API void ParsToK(int Np, int Nk, double *P, double **A, double **Ainv, double *B, double *K);


/*
 And yet another linear transformation, this time to scale initial parameter values to 1.0.
 (In the case where p[i] == 0, we scale by 1 and translate by -1.)
 P = As*S + Bs
 S = As.I * (P - Bs)
*/

extern "C" QUBFAST_API void SetupStartAtOnes(int Np, double *pars, double **Ascal, double *Bscal, double **AscalInv);


#endif