$$\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-2014 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/>.                                        */


#ifdef _WIN32
  #include <windows.h>
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <string.h>
#include <iostream>
#include <fstream>

#define _USE_MATH_DEFINES
#include <math.h>
#include <string.h>
#include <set>
#include <map>
#include <vector>
#include <algorithm>
#include "qubx_model_storage.h"
#include "ublas_matrixutil.h"
#include "qub_constraints.h"
#include "qub_findloops.h"
#include "max_ll_util.h"

using std::cout;
using std::cerr;
using std::endl;
using std::set;
using std::map;
using std::pair;
using std::copy;
using std::min;
using std::max;

qubx_model_storage::qubx_model_storage(int& Nstate_, int& Nclass_, int *class_of_state, callbk_reportfun report_cb_, void *report_data_,
				       int &IeqFv_, double &amp0_, double &var0_, double &Vrev_, double &Vleak_, double &Gleak_ )
  : Nstate(Nstate_), Nclass(Nclass_), IeqFv(IeqFv_), clazz(Nstate_),
    P0(Nstate_), amp(Nclass_), var(Nclass_), d_amp(Nclass_), d_var(Nclass_),
    K0(Nstate_, Nstate_), K1(Nstate_, Nstate_), d_K0(Nstate_, Nstate_), d_K1(Nstate_, Nstate_),
    L(Nstate_, Nstate_), V(Nstate_, Nstate_),
    report_cb(report_cb_), report_data(report_data_), output(report_cb_, report_data_),
    amp0(amp0_), var0(var0_), Vrev(Vrev_), Vleak(Vleak_), Gleak(Gleak_),
    K2(Nstate_, Nstate_), d_K2(Nstate_, Nstate_), P(Nstate_, Nstate_)
{
  if ( class_of_state )
    memcpy(clazz.p, class_of_state, Nstate_*sizeof(int));
}

extern "C" QUBFAST_API qubx_model* qubx_model_create(int Nstate, int Nclass, int* class_of_state,
								      callbk_reportfun report_cb, void *report_data)
{
  qubx_model *m = new qubx_model;
  m->version = 2;
  m->Nstate = Nstate;
  m->Nclass = Nclass;
  qubx_model_storage *st = new qubx_model_storage(m->Nstate, m->Nclass, class_of_state, report_cb, report_data,
						  m->IeqFv, m->amp0, m->var0, m->Vrev, m->Vleak, m->Gleak);
  m->storage = st;
  m->clazz = st->clazz.p;
  m->usePeq = m->IeqFv = 0;
  m->Nchannel = 1;
  m->stop_flag = 0;
  m->progress = 0;
  if ( Nstate > 0 ) {
    m->amp0 = m->var0 = 0.0;
    m->P0 = st->P0.p;
    m->amp = st->amp.p;
    m->var = st->var.p;
    m->d_amp = st->d_amp.p;
    m->d_var = st->d_var.p;
    m->K0 = st->K0.pp;
    m->K1 = st->K1.pp;
    m->K2 = st->K2.pp;
    m->d_K0 = st->d_K0.pp;
    m->d_K1 = st->d_K1.pp;
    m->d_K2 = st->d_K2.pp;
    m->L = st->L.pp;
    m->V = st->V.pp;
    m->P = st->P.pp;
  }
  return m;
}

extern "C" QUBFAST_API void qubx_model_free(qubx_model *m)
{
  delete (qubx_model_storage*) m->storage;
  delete m;
}

#ifdef EMSCRIPTEN
extern "C" EMSCRIPTEN_KEEPALIVE void qubx_model_to_console(qubx_model *m)
{
  EM_ASM_ARGS({
      console.log('Nstate:', $0, 'Nclass:', $1, 'clazz[0]:', $2, 'clazz[1]:', $3, 'P0[0]:', $4, 'P0[1]:', $5);
    }, m->Nstate, m->Nclass, m->clazz[0], m->clazz[1], m->P0[0], m->P0[1]);
  EM_ASM_ARGS({
      console.log('usePeq:', $0, 'IeqFv:', $1, 'Nchannel:', $2, 'iVoltage:', $3);
    }, m->usePeq, m->IeqFv, m->Nchannel, m->iVoltage);
  EM_ASM_ARGS({
      console.log('Gleak:', $0, 'Vleak:', $1, 'Vrev:', $2, 'amp0', $3, 'var0', $4);
    }, m->Gleak, m->Vleak, m->Vrev, m->amp0, m->var0);
  EM_ASM_ARGS({
      console.log('Amp[0]', $0, 'Amp[1]', $1, 'Var[0]', $2, 'Var[1]', $3);
    }, m->amp[0], m->amp[1], m->var[0], m->var[1]);
  EM_ASM_ARGS({
      console.log('k0 00:', $0, 'k0 01:', $1, 'k0 10:', $2, 'k0 11:', $3);
    }, m->K0[0][0], m->K0[0][1], m->K0[1][0], m->K0[1][1]);
  EM_ASM_ARGS({
      console.log('k0 00:', $0, 'k0 01:', $1, 'k0 10:', $2, 'k0 11:', $3);
    }, m->K0[0][0], m->K0[0][1], m->K0[1][0], m->K0[1][1]);
  EM_ASM_ARGS({
      console.log('k1 00:', $0, 'k1 01:', $1, 'k1 10:', $2, 'k1 11:', $3);
    }, m->K1[0][0], m->K1[0][1], m->K1[1][0], m->K1[1][1]);
  EM_ASM_ARGS({
      console.log('l 00:', $0, 'l 01:', $1, 'l 10:', $2, 'l 11:', $3);
    }, m->L[0][0], m->L[0][1], m->L[1][0], m->L[1][1]);
  EM_ASM_ARGS({
      console.log('v 00:', $0, 'v 01:', $1, 'v 10:', $2, 'v 11:', $3);
    }, m->V[0][0], m->V[0][1], m->V[1][0], m->V[1][1]);
  EM_ASM_ARGS({
      console.log('stop:', $0, 'progress:', $1);
    }, m->stop_flag, m->progress);
  EM_ASM_ARGS({
      console.log('&m:', $0, 'L:', $1, 'L[0]:', $2);
    }, m, m->L, m->L[0]);
}
#endif

//--------------------------------------------------------------------------------
// The CChannelCombine class iterates through all combinations of states to make 
// a multi channel model.  Ex( Ch=2, States=3 ) : [2 0 0] [1 1 0] [0 2 0] [1 0 1]
// [0 1 1] [0 0 2].  A c function equivalent is also in qublib nexcom0()
//--------------------------------------------------------------------------------
class CChannelCombine { 
private :
  int m_iChannels;
  int m_iStates;
  int m_iTemp;
  int m_iStateIndex;
  bool m_lContinue;
  int * m_aStates;
public : 
  // Initialize and create the first combination
  CChannelCombine(int iChannels, int iStates);
  // Check if another combination exists 
  bool GetNextCombine();
  // Access the state elements with a 0 based index ( internally stored as 1 based array ) 
  int State0(int i);
  ~CChannelCombine();
};

CChannelCombine :: CChannelCombine(int iChannels, int iStates ) { 
	m_aStates = new int[iStates+1];
	m_iChannels=iChannels;
	m_iStates=iStates;
	m_aStates[1]=iChannels;
	m_iTemp=iChannels;
	m_iStateIndex=0;
	for( int iState=2; iState<=iStates; ++iState)
		m_aStates[iState]=0;
	m_lContinue=true;
	}

bool CChannelCombine :: GetNextCombine( ) { 
	if( m_aStates[m_iStates] == m_iChannels )
		return false;
	if(m_iTemp > 1) 
		m_iStateIndex=0;
	m_iStateIndex++;
	m_iTemp = m_aStates[m_iStateIndex];
	m_aStates[m_iStateIndex] = 0;
	m_aStates[1] = m_iTemp - 1;
	m_aStates[m_iStateIndex + 1] = m_aStates[m_iStateIndex + 1] + 1;
	return true;
	}

int CChannelCombine :: State0(int i){
	return m_aStates[i+1];
	}

CChannelCombine :: ~CChannelCombine() { 
	delete [] m_aStates;
	}


//-----------------------------------------------------------------------
// Slightly Similar to a numerical recipes algorithm eclazz, but this 
// returns a sorted aClasses array (though it is less efficient).
// aData - array of values.  iData - number of values.
// iClasses - output number of classes.  aClasses - output class values.
// ia[0..iData-1] - output  index to aClasses for each aData[0..iData-1]
// This function iteratively finds the minimum value in the list, then 
// considers all values within 1.0e-5 of that to be equivalent class and 
// removes them from the list until all data have been placed in classes.
void eclass(float * aData, int iData, int & iClasses, float * aClasses, 
			int * ia, double openDirection) {
	int  i,k;
	float amin;

	int iCopyPos;				// While scanning aClasses array, position being copied to ( = current - # of elements at current class already processed )
	bool lCompare = ( openDirection > 0 );	// whether to find min or max ( ascending or descending )

	int * aRemain = new int[iData];	// Array of indices of items to be classified
	for (i=0; i<iData; i++) 
		aRemain[i] = i;

	iClasses = 0;
	while( iData > 0 ) {
		// Find the minimum (or maximum ) value in the remaining aRemain list 
		k = aRemain[0];
		amin = aData[k];
		
		for( i=1; i<iData; i++ ) {
			k=aRemain[i];
			if( lCompare == (aData[k]<amin)) 
				amin=aData[k];
			}
		
		// add amin to aClasses[] array of classes and all aData[k]~= amin set ia[k]=iClasses
		aClasses[iClasses]=amin;
		iCopyPos=0;
		for (i=0; i<iData; i++) {
			k = aRemain[i]; 
			if (fabs(aData[k]-amin)<1.0e-5) 
				ia[k]=iClasses;					// include as result for this minimum ...
			else 
				aRemain[iCopyPos++]=k;				// ... or copy back to array for next iteration
			}
		iData=iCopyPos;

		iClasses++;
		}   
	
	delete [] aRemain;
	}



//-----------------------------------------------------------------------------
// multi-channel Markov model creation : this function will make a multichannel 
// markov model out of a single channel markov model
// Input nchannel, nstate, npath, state, path, ratio
// Output nmutistate, nmutipath, nmuticlass, mutistate, mutipath
//-----------------------------------------------------------------------
extern "C" QUBFAST_API void qubx_mutimakv_size(int nchannel, int nstate, QUBFAST_VAR_NOT_USED int npath, QUBFAST_VAR_NOT_USED int *path,
					       int *nmutistate, int *nmutipath)
// returns upper bound on array sizes needed for mutistate and mutipath, via args
{
	int i, iPaths=0, iState;
	
	//-- Precalc model size for memory allocation 
	// Calculates the number of Combinations of channels/states per row in mutistate.
	// For each state, add to temp log(state+nchann)-log(state) and then antilog on the sum
	//float tmp=0.0;
	//for( n=1; n<nstate; n++) 
	  //	tmp += float(log(double(n+nchannel))-log(double(n))); 
	//n = int(exp(tmp)+1);						//estimates nMutiState - may be 1 greater than necc.
	//*nmutistate = n * (nstate+1);


	//----- Fill in mutistate , acond, iPaths.
	// mutistate is an integer array of (estimated nMutiState) * (nstate+1)
	// Each (nState+1) group is a 1-based array of number of channels in each state,
	// and [0] element is set to the eclass result class for the row
	// iterate through combinations of states / channels and set nMutiState
	// aCond is the effective conductance <aCond[]> = sum[all states] ( channels * ratio ) 
	// iPaths is the number of mutistate connections which have only one channel transition
	*nmutistate = 0;
	CChannelCombine oMeta(nchannel, nstate);
	do {
		for ( i=0; i<nstate; ++i ) {
			iState=oMeta.State0(i);
			if( iState > 0 )								// pre calc # of paths 
				iPaths += (nstate-1);
			}
		++(*nmutistate);
		} while( oMeta.GetNextCombine() );

	//-----
	// precalc the number of transition paths for memory allocation
	// iterates through all combinations of mutistates and counts those 
	// which differ only by 1 transition.
	// New method : for each non zero state, add paths for the number of 
	//  places one of the channels could go.  This is precalced above.
	// This calculates a maximum possible, but the actual number used is less
	/*
	n=0;
	for(k=0; k<nmutistate; k++)  
		for (l=0; l<nmutistate; l++) { 
			s=0;
			for (i=0; i<nstate; i++)  
				s += int(fabs(aMutiState[k][i+1]-aMutiState[l][i+1])); 
			if (s==2) 
				n++; 
			} 
	mutipath = new int[n*4];
	*/
	*nmutipath = iPaths;
}

extern "C" QUBFAST_API void qubx_mutimakv(int nchannel, int nstate, int npath, int *state, int *path, float *ratio, 
					  int nmutistate, int nmutipath, int* nmuticlass, int* mutistate, int* mutipath)
// mutistate = new int[nmutistate * (nstate + 1)]
// mutipath = new int[nmutipath * 4]
{
	int i,k=0,l=0,m,n,s,t, iPaths=0, iState;
	
    //-- Determine whether open is negative or positive and use to set sort order 
	// for muticlasses in eclass().
	// Note : The try catch is probably not effective here.  Accessing beyond the 
	// bounds of the array would likely run into more data and not throw an exception.
	// Better to ensure that the passed array ratio[] has an element [1] always.
	double openDirection;
	try {			
		openDirection = ratio[1];
		}
	catch (...) {
		openDirection = 1.0;
		}

	float * aCond = new float[nmutistate];		// Conductance class sums
	
	// Create a 2d style array for accessing mutistate array 
	int ** aMutiState = new int * [nmutistate];
	for( i=0; i<nmutistate; ++i ) 
		aMutiState[i] = &(mutistate[i*(nstate+1)]);

	//----- Fill in mutistate , acond, iPaths.
	// mutistate is an integer array of (estimated nMutiState) * (nstate+1)
	// Each (nState+1) group is a 1-based array of number of channels in each state,
	// and [0] element is set to the eclass result class for the row
	// iterate through combinations of states / channels and set nMutiState
	// aCond is the effective conductance <aCond[]> = sum[all states] ( channels * ratio ) 
	// iPaths is the number of mutistate connections which have only one channel transition
	nmutistate=0;
	CChannelCombine oMeta(nchannel, nstate);
	do {
		aCond[nmutistate]=0.0;
		for ( i=0; i<nstate; ++i ) {
			iState=oMeta.State0(i);
			aMutiState[nmutistate][i+1]=iState;
			aCond[nmutistate] += ((float) iState) * ratio[state[i]];
			if( iState > 0 )								// pre calc # of paths 
				iPaths += (nstate-1);
			}
		++nmutistate;
		} while( oMeta.GetNextCombine() );

	//----- conductance classes 
	float * aClasses = new float[nmutistate];	// output buffer, discarded unused
	int *  ia = new int[nmutistate];			// mapping of data to classes
	
	eclass( aCond, nmutistate, *nmuticlass, aClasses, ia, openDirection);		// Fill C,ia with classes
	
	for (n=0; n<nmutistate; n++) 
		aMutiState[n][0] = ia[n];
	
	delete [] aCond;
	delete [] aClasses;
	delete [] ia;

	// transition paths - for each mutipath pair (m,n), if they differ 
	// by only 1 channel state change then find the original states (k,l)
	// and check if a path exists between k  & l -->> add a mutipath if so.
	//
	// Note : Consider a transition where  two channels change state in a 
	// single delta T - is this accounted for elsewhere ? {JB}
	// In a simple 2 state 3 channel model currently (C,C,C) we are only 
	// setting mutipaths for (C,C,O),(C,O,C),(O,C,C) but all states 
	// including (O,O,O) are reachable in a single delta T ?
	//
	// Ans: A=exp(Q*dt) doesn't have any zeros -- it allows for many quick transitions within one dt. {CN}
	nmutipath = 0; 
	for( m=0; m<nmutistate; m++ ) {
		for( n=0; n<nmutistate; n++ ) {
			s=0;
			for( i=0; i<nstate; i++) 
				s += abs(aMutiState[m][i+1]-aMutiState[n][i+1]);
			if (s==2) {	// if muti[m] can become muti[n] by only 1 channel changing state 
				// Find k,l the source and dest states 
				for (i=0; i<nstate; i++) {
					t = aMutiState[m][i+1]-aMutiState[n][i+1];
					if (t==1) 
						k = i;
					if (t==-1) 
						l = i;
					}
				// Check for a path from k to l in model 
				for (i=0; i<npath; i++) {
					if (path[i*4+0]==k && path[i*4+1]==l) {   
						// add a nutipath 
						mutipath[nmutipath*4+0] = m;
						mutipath[nmutipath*4+1] = n;
						mutipath[nmutipath*4+2] = k;
						mutipath[nmutipath*4+3] = l;
						nmutipath++;
						}
					}
				}
			}
		}

	delete [] aMutiState;
	}


void ktomutik(int nstate, QUBFAST_VAR_NOT_USED int npath, QUBFAST_VAR_NOT_USED int* path, int nmutistate, 
	      int* mutistate, int nmutipath, int* mutipath, double *k0, double *k1, double *k2, int *L, int *V, int *P,
	      double *mutik0, double *mutik1, double *mutik2, int *mutiL, int *mutiV, int *mutiP) {
	int  I,J,i,j,m;
    
	// rate constants 
	for (I=0; I<nmutistate; I++) {
		for (J=0; J<nmutistate; J++) { 
			mutik0[I*nmutistate+J] = 0.0;
			mutik1[I*nmutistate+J] = 0.0;
			mutik2[I*nmutistate+J] = 0.0;
			mutiL[I*nmutistate+J] = 0;
			mutiV[I*nmutistate+J] = 0;
			mutiP[I*nmutistate+J] = 0;
		}
	}

	for (m=0; m<nmutipath; m++) {
		I = mutipath[m*4+0];
		J = mutipath[m*4+1];
		i = mutipath[m*4+2];
		j = mutipath[m*4+3];
		mutik0[I*nmutistate+J] = mutistate[I*(nstate+1)+i+1]*k0[i*nstate+j];
		mutik1[I*nmutistate+J] = k1[i*nstate+j];
		mutik2[I*nmutistate+J] = k2[i*nstate+j];
		mutiL[I*nmutistate+J] = L[i*nstate+j];
		mutiV[I*nmutistate+J] = V[i*nstate+j];
		mutiP[I*nmutistate+J] = P[i*nstate+j];
		}
}


qubx_multi_model_storage::qubx_multi_model_storage(int Nstate, int& npath_, int& nmutistate_, int& nmutipath_, int& nmuticlass_)
  : nmutistate(nmutistate_), nmutipath(nmutipath_), nmuticlass(nmuticlass_), npath(npath_), class_order(Nstate)
{}

void qubx_multi_model_storage::setup(int nmutistate_, int nmutipath_, qubx_multi_model *mm)
{
  nmutistate = nmutistate_;
  nmutipath = nmutipath_;
  cmutistate.resize(nmutistate_);
  mutistate.resize(nmutistate_, mm->single_model->Nstate+1);
  mutipath.resize(nmutipath_, 4);
  mm->cmutistate = cmutistate.p;
  mm->mutistate = mutistate.p;
  mm->mutipath = mutipath.p;
}
  
extern "C" QUBFAST_API qubx_multi_model* qubx_multi_model_create(qubx_model *single_model)
{
  qubx_multi_model *mm = new qubx_multi_model;
  mm->version = 1;
  mm->single_model = single_model;
  qubx_multi_model_storage *st = new qubx_multi_model_storage(single_model->Nstate, mm->npath, mm->nmutistate, mm->nmutipath, mm->nmuticlass);
  mm->storage = st;
  mm->class_order = st->class_order.p;

  int nstate = single_model->Nstate;
  
  // renumber classes to be contiguous starting at 0:
  std::vector<int> occupied(single_model->Nclass);
  for (int i=0; i<nstate; ++i)
    occupied[single_model->clazz[i]] = 1;
  int nclass = 0;
  for (int i=0; i<single_model->Nclass; ++i) {
    if ( occupied[i] ) {
      mm->class_order[nclass] = i;
      occupied[i] = nclass++;
    }
  }
  mm->nclass_order = nclass;
  std::vector<int> clazz(nstate);
  for (int i=0; i<nstate; ++i)
    clazz[i] = occupied[single_model->clazz[i]];

  // derive legacy path, npath
  int npath = 0;
  for (int a=0; a<nstate; ++a) {
    for (int b=a+1; b<nstate; ++b) {
      if ( single_model->K0[a][b] )
	npath += 2;
    }
  }
  st->npath = npath;
  st->path.resize(npath, 4);
  mm->path = st->path.p;
  int ipath = 0;
  for (int a=0; a<nstate; ++a) {
    for (int b=a+1; b<nstate; ++b) {
      if ( single_model->K0[a][b] ) {
	st->path.m(ipath, 0) = st->path.m(ipath+1, 1) = a;
	st->path.m(ipath, 1) = st->path.m(ipath+1, 0) = b;
	ipath += 2;
      }
    }
  }


  // ratio: excess class amplitude (minus baseline)
  fq::vector<float> ratio( nclass );
  float amp0 = (float) (single_model->iVoltage ? 0.0 : single_model->amp[0]);
  for (int i=0; i<nclass; ++i )
    ratio[i] = ((float) single_model->amp[mm->class_order[i]]) - amp0;
  
  qubx_mutimakv_size(single_model->Nchannel, nstate, npath, mm->path, &mm->nmutistate, &mm->nmutipath);
  st->setup(mm->nmutistate, mm->nmutipath, mm);
  
  qubx_mutimakv( single_model->Nchannel, nstate, npath, &clazz[0], mm->path, ratio, 
		 mm->nmutistate, mm->nmutipath, &mm->nmuticlass, mm->mutistate, mm->mutipath);

  for (int i=0; i<mm->nmutistate; ++i )
    mm->cmutistate[i] = mm->mutistate[i*(nstate+1)+0]; // mutistate[i][0] -- muticlass of mutistate i
  
  qubx_model_storage *single_st = (qubx_model_storage *) single_model->storage;
  qubx_model *multi_model = qubx_model_create(mm->nmutistate, mm->nmuticlass, mm->cmutistate, single_st->report_cb, single_st->report_data);
  mm->multi_model = multi_model;
  multi_model->usePeq = single_model->usePeq;
  multi_model->Nchannel = 1;
  multi_model->Gleak = single_model->Gleak;
  multi_model->Vleak = single_model->Vleak;
  multi_model->Vrev = single_model->Vrev;
  memcpy(multi_model->clazz, mm->cmutistate, mm->nmutistate*sizeof(int));

  return mm;
}

extern "C" QUBFAST_API void qubx_multi_model_free(qubx_multi_model *mm)
{
  qubx_model_free(mm->multi_model);
  delete (qubx_multi_model_storage*) mm->storage;
  delete mm;
}

extern "C" QUBFAST_API void qubx_multi_model_update(qubx_multi_model *mm)
{
	int nstate = mm->single_model->Nstate;

        ktomutik(nstate, mm->npath, mm->path, mm->nmutistate, mm->mutistate, mm->nmutipath, 
		 mm->mutipath, mm->single_model->K0[0], mm->single_model->K1[0], mm->single_model->K2[0],
		 mm->single_model->L[0], mm->single_model->V[0], mm->single_model->P[0],
		 mm->multi_model->K0[0], mm->multi_model->K1[0], mm->multi_model->K2[0],
		 mm->multi_model->L[0], mm->multi_model->V[0], mm->multi_model->P[0]);

	// state entry probability is not dealt with here...

	int Nclass = mm->nclass_order;
	fq::vector<int> clazz(nstate);
	fq::vector<int> cls_renumber(mm->single_model->Nclass);
	for (int i=0; i<Nclass; ++i)
	  cls_renumber[ mm->class_order[i] ] = i;
	for (int i=0; i<nstate; ++i)
	  clazz[i] = cls_renumber[ mm->single_model->clazz[i] ];

	fq::vector<double> exc_amp(Nclass);
	fq::vector<double> exc_var(Nclass);
	double amp0, var0;
	if ( mm->single_model->iVoltage ) {
	  mm->multi_model->IeqFv = mm->single_model->IeqFv;
	  mm->multi_model->iVoltage = mm->single_model->iVoltage;
	  mm->multi_model->amp0 = mm->single_model->amp0;
	  mm->multi_model->var0 = mm->single_model->var0;
	  amp0 = var0 = 0.0;
	  for (int i=0; i<Nclass; ++i) {
	    exc_amp[i] = mm->single_model->amp[ mm->class_order[i] ];
	    exc_var[i] = mm->single_model->var[ mm->class_order[i] ];
	  }
	} else {
	  amp0 = mm->multi_model->amp0 = mm->single_model->amp[0];
	  var0 = mm->multi_model->var0 = mm->single_model->var[0];
	  for (int i=0; i<Nclass; ++i) {
	    exc_amp[i] = mm->single_model->amp[ mm->class_order[i] ] - amp0;
	    exc_var[i] = mm->single_model->var[ mm->class_order[i] ] - var0;
	  }
	}

	double *mutiamp = mm->multi_model->amp;
	double *mutivar = mm->multi_model->var;

	for (int i=0; i<mm->nmuticlass; i++) {
		mutiamp[i] = amp0; // closed_amp + sum(open_amp === ratio)
		//mutixms[i] = 0.0;          // sqrt( sum( xms*xms ) )
		// now this should probably be // sqrt( closed_noise^2 + sum(excess_noise^2) ) ??
		mutivar[i] = var0; // will be sqrt( closed_noise^2 [var0] + sum(excessNoise^2) )
		//float sum=0.0; // # states contributing to this mutistate; should always be == nchannel

                int k;
		for (k=0; mm->mutistate[k*(nstate+1)+0]!=i; k++)    // mutistate[k][0]
			; // find the first mutistate in this muticlass.
		      // if (e.g. 2 sub-opens == 1 open 1 closed)
		      //   then the first example may not have representative s.d. and allXFixed
		
		for (int j=0; j<nstate; j++) {                      // mutistate[k][j+1]
			int countStateJ = mm->mutistate[k*(nstate+1)+j+1];
			mutiamp[i] += countStateJ*exc_amp[clazz[j]];
			mutivar[i] += countStateJ*exc_var[clazz[j]];
		}
	}

}

qubx_constrained_storage::qubx_constrained_storage(int& Npar_, int& NparR_, int& NfPar_, int& Ncns_, int& NcnsR_)
  : Npar(Npar_), NparR(NparR_), NfPar(NfPar_), Ncns(Ncns_), NcnsR(NcnsR_),
    Ain(3*Npar, Npar), Acns(Npar, Npar), Ainv(Npar, Npar), Ascal(Npar, Npar), Asci(Npar, Npar),
    Bin(3*Npar), Bcns(Npar), Bscal(Npar), pars(Npar), fPars_unscaled(Npar), fPars(Npar),
    kR_indices(Npar_)
{
  Ncns = NfPar = 0;
}

int qubx_constrained_storage::setup()
{
  if ( ! NfPar ) {
    constraint_indices.resize(Ncns);
    IdentifyFixedParams(Npar, Ncns, pars.p, Ain.pp, Bin.p, &NparR, kR_indices.p, &NcnsR, constraint_indices.p);
    parsR.resize(NparR);
    AinR.resize(NcnsR, NparR);
    BinR.resize(NcnsR);
    EliminateFixedParams(Npar, Ncns, pars.p, Ain.pp, Bin.p, NparR, kR_indices.p, NcnsR, constraint_indices.p, parsR.p, AinR.pp, BinR.p);
    ReduceConstraints(NparR, &NcnsR, AinR.pp, BinR.p);
    NfPar = NparR - NcnsR;
  }
  if ( NcnsR < NparR ) {
    int rtn = SetupLinearConstraints(NcnsR, AinR.pp, BinR.p, NparR, parsR.p, fPars_unscaled.p,
				     Acns.pp, Bcns.p, Ainv.pp, NULL, NULL);
    if ( 0 == rtn ) {
      SetupStartAtOnes(NfPar, fPars_unscaled.p, Ascal.pp, Bscal.p, Asci.pp);
      for (int i=0; i<NfPar; ++i)
	fPars.p[i] = 1.0;
    }
    return rtn;
  }
  return NfPar;
}
  
void qubx_constrained_storage::free_to_pars()
{
  if ( NcnsR < NparR ) {
    ParsToK(NfPar, NfPar, fPars.p, Ascal.pp, Asci.pp, Bscal.p, fPars_unscaled.p);
    ParsToK(NfPar, NparR, fPars_unscaled.p, Acns.pp, Ainv.pp, Bcns.p, parsR.p);
    Kr_to_K(parsR.p, NparR, pars.p, kR_indices.p);
  }
}
  
void qubx_constrained_storage::pars_to_free()
{
  if ( Ncns < Npar ) {
    for (int i=0; i<NparR; ++i)
      parsR.p[i] = pars.p[kR_indices.p[i]];
    KtoPars(NparR, NfPar, parsR.p, Acns.pp, Ainv.pp, Bcns.p, fPars_unscaled.p);
    KtoPars(NfPar, NfPar, fPars_unscaled.p, Ascal.pp, Asci.pp, Bscal.p, fPars.p);
  }
}

void qubx_constrained_storage::calc_std(double *InvHessian, double *dPar)
{  
  matrix<double> Acombined = prod(Acns.m, Ascal.m);
  int n = NparR;
  int m = NfPar;
  // each element squared:
  for (int i=0; i<n; ++i)
    for (int j=0; j<m; ++j)
      Acombined(i,j) *= Acombined(i,j);
  // mil confidence: 
  //   dx[i] = sqrt( sum_j( A[i,j]**2 * H[j,j] ) )    { * exp(x[i]) }
  for (int i=0; i<n; ++i) {
    double sum = 0.0;
    for (int j=0; j<m; ++j) {
      sum += Acombined(i,j) * InvHessian[j*m+j];
    }
    dPar[kR_indices.p[i]] = sum;
  }
}
   

extern "C" QUBFAST_API qubx_constrained* qubx_constrained_create(int Npar)
{
  qubx_constrained *obj = new qubx_constrained;
  obj->Npar = Npar;
  qubx_constrained_storage *st = new qubx_constrained_storage(obj->Npar, obj->NparR, obj->NfPar, obj->Ncns, obj->NcnsR);
  obj->version = 1;
  obj->storage = st;
  obj->Ain = st->Ain.pp;
  obj->Bin = st->Bin.p;
  obj->pars = st->pars.p;
  obj->fPars = st->fPars.p;
  return obj;
}

extern "C" QUBFAST_API void qubx_constrained_free(qubx_constrained *obj)
{
  qubx_constrained_storage *st = (qubx_constrained_storage *) obj->storage;
  delete st;
  delete obj;
}

extern "C" QUBFAST_API int qubx_constrained_setup(qubx_constrained *obj)
{
  qubx_constrained_storage *st = (qubx_constrained_storage *) obj->storage;
  int rtn = st->setup();
  obj->NfPar = st->NfPar;
  obj->Ncns = st->Ncns;
  obj->Acns = st->Acns.pp;
  obj->Ascal = st->Ascal.pp;
  obj->Bcns = st->Bcns.p;
  obj->Bscal = st->Bscal.p;
  return rtn;
}

extern "C" QUBFAST_API void qubx_constrained_free_to_pars(qubx_constrained *obj)
{
  qubx_constrained_storage *st = (qubx_constrained_storage *) obj->storage;
  st->free_to_pars();
}

extern "C" QUBFAST_API void qubx_constrained_pars_to_free(qubx_constrained *obj)
{
  qubx_constrained_storage *st = (qubx_constrained_storage *) obj->storage;
  st->pars_to_free();
}

extern "C" QUBFAST_API void qubx_constrained_calc_std(qubx_constrained *obj, double *InvHessian, double *dPar)
{
  qubx_constrained_storage *st = (qubx_constrained_storage *) obj->storage;
  st->calc_std(InvHessian, dPar);
}

    
qubx_stim_amps_storage::qubx_stim_amps_storage(qubx_model *m, int Nstim_, int Nsig_, double *stimclasses_, double *sc_frac_)
  : model(m), Nstim(Nstim_), Nsig(Nsig_),
    stimclasses(Nstim_*Nsig_), sc_frac(Nstim_),
    sc_amp(Nstim_, m->Nclass), sc_var(Nstim_, m->Nclass), sc_std(Nstim_, m->Nclass)
{
  memcpy(stimclasses.p, stimclasses_, Nstim*Nsig*sizeof(double));
  memcpy(sc_frac.p, sc_frac_, Nstim*sizeof(double));
  setup_scs();
  
  cns = qubx_constrained_create((m->IeqFv ? 2 : 0) + 2*m->Nclass);
  cnst = (qubx_constrained_storage *) cns->storage;
  model_to_pars();
}

qubx_stim_amps_storage::~qubx_stim_amps_storage()
{
  qubx_constrained_free(cns);
}
  
void qubx_stim_amps_storage::setup_scs()
{
  int iVoltage = model->IeqFv ? model->iVoltage : 0;
  if ( iVoltage ) {
    for (int sc=0; sc<Nstim; ++sc) {
      double volts = 1e-3 * stimclasses.p[sc * Nsig + iVoltage];
      double dV = volts - 1e-3 * model->Vrev;
      double Ileak = (volts - 1e-3*model->Vleak)*model->Gleak;
      for (int c=0; c<model->Nclass; ++c) {
	sc_amp.m(sc,c) = model->amp0 + dV*model->amp[c] + Ileak;
	sc_var.m(sc,c) = model->var0 + dV*dV*model->var[c];
	sc_std.m(sc,c) = sqrt(sc_var.m(sc,c));
      }
    }
  } else {
    for (int sc=0; sc<Nstim; ++sc) {
      memcpy(sc_amp.pp[sc], model->amp, model->Nclass*sizeof(double));
      memcpy(sc_var.pp[sc], model->var, model->Nclass*sizeof(double));
      for (int c=0; c<model->Nclass; ++c)
	sc_std.m(sc,c) = (model->var[c] > 0.0) ? sqrt(model->var[c]) : 0.0;
    }
  }
}

void qubx_stim_amps_storage::pars_to_model()
{
  int p=0;
  for (int i=0; i<model->Nclass; ++i) {
    model->amp[i] = cns->pars[p++];
    model->var[i] = fabs(cns->pars[p++]);
  }
  if ( model->IeqFv && model->iVoltage ) {
    model->amp0 = cns->pars[p++];
    model->var0 = fabs(cns->pars[p++]);
  }
}

void qubx_stim_amps_storage::model_to_pars()
{
  int p=0;
  for (int i=0; i<model->Nclass; ++i) {
    cns->pars[p++] = model->amp[i];
    cns->pars[p++] = model->var[i];
  }
  if ( model->IeqFv && model->iVoltage ) {
    cns->pars[p++] = model->amp0;
    cns->pars[p++] = model->var0;
  }
}

extern "C" QUBFAST_API qubx_stim_amps*  qubx_stim_amps_create(qubx_model *m, int Nstim_, int Nsig_,
							     double *stimclasses_, double *sc_frac_)
{
  qubx_stim_amps *obj = new qubx_stim_amps;
  qubx_stim_amps_storage *st = new qubx_stim_amps_storage(m, Nstim_, Nsig_, stimclasses_, sc_frac_);
  obj->version = 1;
  obj->storage = st;
  obj->model = m;
  obj->Nstim = st->Nstim;
  obj->Nsig = st->Nsig;
  obj->stimclasses = st->stimclasses.p;
  obj->sc_frac = st->sc_frac.p;
  obj->cns = st->cns;
  obj->sc_amp = st->sc_amp.pp;
  obj->sc_std = st->sc_std.pp;
  return obj;
}

extern "C" QUBFAST_API void qubx_stim_amps_free(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  delete st;
  delete obj;
}
  
extern "C" QUBFAST_API void qubx_stim_amps_setup_scs(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  st->setup_scs();
}

extern "C" QUBFAST_API void qubx_stim_amps_pars_to_model(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  st->pars_to_model();
}

extern "C" QUBFAST_API void qubx_stim_amps_model_to_pars(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  st->model_to_pars();
}

extern "C" QUBFAST_API void qubx_stim_amps_fix_amp0(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  qubx_constrained *c = obj->cns;
  for (int j=0; j<c->Npar; ++j)
    c->Ain[c->Ncns][j] = 0.0;
  c->Ain[c->Ncns][2*st->model->Nclass] = 1.0;
  c->Bin[c->Ncns] = st->model->amp0;
  ++c->Ncns;
}

extern "C" QUBFAST_API void qubx_stim_amps_fix_var0(qubx_stim_amps *obj)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  qubx_constrained *c = obj->cns;
  for (int j=0; j<c->Npar; ++j)
    c->Ain[c->Ncns][j] = 0.0;
  c->Ain[c->Ncns][2*st->model->Nclass+1] = 1.0;
  c->Bin[c->Ncns] = st->model->var0;
  ++c->Ncns;
}

extern "C" QUBFAST_API void qubx_stim_amps_fix_amp(qubx_stim_amps *obj, int cls)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  qubx_constrained *c = obj->cns;
  for (int j=0; j<c->Npar; ++j)
    c->Ain[c->Ncns][j] = 0.0;
  c->Ain[c->Ncns][2*cls] = 1.0;
  c->Bin[c->Ncns] = st->model->amp[cls];
  ++c->Ncns;
}

extern "C" QUBFAST_API void qubx_stim_amps_fix_var(qubx_stim_amps *obj, int cls)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  qubx_constrained *c = obj->cns;
  for (int j=0; j<c->Npar; ++j)
    c->Ain[c->Ncns][j] = 0.0;
  c->Ain[c->Ncns][2*cls+1] = 1.0;
  c->Bin[c->Ncns] = st->model->var[cls];
  ++c->Ncns;
}

extern "C" QUBFAST_API void qubx_stim_amps_calc_std(qubx_stim_amps *obj, double *InvHessian)
{
  qubx_stim_amps_storage *st = (qubx_stim_amps_storage *) obj->storage;
  std::vector<double> dPar(obj->cns->Npar);
  qubx_constrained_calc_std(obj->cns, InvHessian, & dPar[0]);

  int p=0;
  for (int i=0; i<st->model->Nclass; ++i) {
    st->model->d_amp[i] = dPar[p++];
    st->model->d_var[i] = dPar[p++];
  }
  if ( st->model->IeqFv && st->model->iVoltage ) {
    st->model->d_amp0 = dPar[p++];
    st->model->d_var0 = dPar[p++];
  }
}

qubx_stim_rates_storage::qubx_stim_rates_storage(qubx_model *m, int Nstim_, int Nsig_, double *stimclasses_, double *sc_frac_,
						 double sampling_, eigen_func custom_eigen_)
  : model(m), custom_eigen(custom_eigen_), Nstim(Nstim_), Nsig(Nsig_),
    stimclasses(Nstim_*Nsig_), sc_frac(Nstim), sampling(sampling_),
    QQ(Nstim_),   AA(Nstim_),   QQ_evec(Nstim_),   QQ_einv(Nstim_),   QQ_eig(Nstim_),   QQ_imag(Nstim_),
    QQ_p(Nstim_), AA_p(Nstim_), QQ_evec_p(Nstim_), QQ_einv_p(Nstim_), QQ_eig_p(Nstim_), QQ_imag_p(Nstim_),
    A(m->Nstate, m->Nstate), tmp_diag(m->Nstate, m->Nstate), tmp_m(m->Nstate, m->Nstate)
{
  memcpy(stimclasses.p, stimclasses_, Nstim*Nsig*sizeof(double));
  if ( sc_frac_ )
    memcpy(sc_frac.p, sc_frac_, Nstim*sizeof(double));
  cns = qubx_constrained_create(CountK(m->Nstate));
  cnst = (qubx_constrained_storage *) cns->storage;
  int Ns = model->Nstate;
  for (int sc=0; sc<Nstim; ++sc) {
    QQ[sc].resize(Ns,Ns);          QQ_p[sc]       = QQ[sc].pp;
    AA[sc].resize(Ns,Ns);          AA_p[sc]       = AA[sc].pp;
    QQ_evec[sc].resize(Ns,Ns);     QQ_evec_p[sc]  = QQ_evec[sc].pp;
    QQ_einv[sc].resize(Ns,Ns);     QQ_einv_p[sc]  = QQ_einv[sc].pp;
    QQ_eig[sc].resize(Ns);         QQ_eig_p[sc]   = QQ_eig[sc].p;
    QQ_imag[sc].resize(Ns);        QQ_imag_p[sc]  = QQ_imag[sc].p;
  }
  setup_QQ();
  setup_AA();
}

qubx_stim_rates_storage::~qubx_stim_rates_storage()
{
  qubx_constrained_free(cns);
}

void qubx_stim_rates_storage::setup_QQ()
{
  qubx_model_storage *st = (qubx_model_storage *) model->storage;
  for (int sc=0; sc<Nstim; ++sc)
    BuildQ_p(model->Nstate, QQ[sc].m, st->K0.m, st->K1.m, st->K2.m, st->L.m, st->V.m, st->P.m, stimclasses.p + sc*Nsig);
}

void qubx_stim_rates_storage::setup_AA()
{
  qubx_model_storage *st = (qubx_model_storage *) model->storage;
  bool err_eigen = false;
  if ( custom_eigen ) {
    for (int sc=0; sc<Nstim; ++sc) {
      err_eigen = err_eigen || custom_eigen(model->Nstate, QQ[sc].p, QQ_eig[sc].p, QQ_imag[sc].p, QQ_evec[sc].p, QQ_einv[sc].p);
      if ( sampling )
	expm_eig(QQ_eig[sc].p, QQ_evec[sc].m, QQ_einv[sc].m, sampling, tmp_diag.m, tmp_m.m, AA[sc].m);
    }
  } else {
    int n = model->Nstate;

#pragma omp parallel
    {
      matrix<double> tmpd(n,n), tmpm(n,n);
      tmpd.clear();
#pragma omp for
      for (int sc=0; sc<Nstim; ++sc) {
	eigen_vals_vecs(QQ[sc].m, QQ_eig[sc].p, QQ_imag[sc].p, QQ_evec[sc].m, QQ_einv[sc].m, st->output);
	if ( sampling )
	  expm_eig(QQ_eig[sc].p, QQ_evec[sc].m, QQ_einv[sc].m, sampling, tmpd, tmpm, AA[sc].m);
      }
    }

  }
  if ( err_eigen ) {
    st->output << "No convergence in eigen calculations" << endl;
  }
}

void qubx_stim_rates_storage::setup_A(int sc, double dt)
{
  setup_A(sc, dt, tmp_diag.m, tmp_m.m, A.m);
}

void qubx_stim_rates_storage::setup_A(int sc, double dt,
				      matrix<double>& tmp_diag_x, matrix<double>& tmp_m_x, matrix<double>& A_result)
{
  expm_eig(QQ_eig[sc].p, QQ_evec[sc].m, QQ_einv[sc].m, dt, tmp_diag_x, tmp_m_x, A_result);
}

void qubx_stim_rates_storage::fixRate(int a, int b)
{
  AddFixRateConstraint(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin, a, b);
}

void qubx_stim_rates_storage::fixExp(int a, int b)
{
  AddFixExpConstraint(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin, a, b);
}

void qubx_stim_rates_storage::scaleRate(int a, int b, int c, int d)
{
  AddScaleRateConstraint(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin, a, b, c, d);
}

void qubx_stim_rates_storage::scaleExp(int a, int b, int c, int d)
{
  AddScaleExpConstraint(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin, a, b, c, d);
}

void qubx_stim_rates_storage::balanceLoop(int Nloop, int *loopstates)
{
  AddLoopBalConstraint(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin, Nloop, loopstates, consoleLog, 0);
}

void qubx_stim_rates_storage::balanceLoops()
{
  int Nstate = model->Nstate;
  int Nrate = 0;
  std::vector<int> edges;
  for (int a=0; a<Nstate; ++a) {
    for (int b=a+1; b<Nstate; ++b) {
      if ( model->K0[a][b] ) {
        edges.push_back(a);
        edges.push_back(b);
        edges.push_back(b);
        edges.push_back(a);
        Nrate += 2;
      }
    }
  }
  std::vector<int> loopsizes(Nrate);
  std::vector<int> loopstates(Nrate*Nrate);
  int Nloop = qub_findloops(Nstate, Nrate, & edges[0], & loopsizes[0], & loopstates[0]);
  for (int i=0; i<Nloop; ++i)
    balanceLoop(loopsizes[i], & loopstates[i*Nrate]);
}

int qubx_stim_rates_storage::setup_constraints()
{
  if ( ! cns->NfPar )
    AddAutoConstraints(model->Nstate, model->K0, model->K1, model->K2, model->L, model->V, model->P, & cns->Ncns, cns->Ain, cns->Bin);
  model_to_pars();
  return qubx_constrained_setup(cns);
}

void qubx_stim_rates_storage::pars_to_model()
{
  KtoK012(model->Nstate, cns->pars, model->K0, model->K1, model->K2);
  setup_QQ();
  setup_AA();
}

void qubx_stim_rates_storage::model_to_pars()
{
  K012toK(model->Nstate, model->K0, model->K1, model->K2, cns->pars);
}

extern "C" QUBFAST_API qubx_stim_rates *qubx_stim_rates_create(qubx_model *model, int Nstim_, int Nsig_,
						   double *stimclasses_, double *sc_frac_, double sampling,
						   eigen_func custom_eigen)
{
  qubx_stim_rates *obj = new qubx_stim_rates;
  qubx_stim_rates_storage *st = new qubx_stim_rates_storage(model, Nstim_, Nsig_, stimclasses_, sc_frac_,
							    sampling, custom_eigen);
  obj->version = 1;
  obj->storage = st;
  obj->model = model;
  obj->custom_eigen = custom_eigen;
  obj->Nstim = st->Nstim;
  obj->Nsig = st->Nsig;
  obj->stimclasses = st->stimclasses.p;
  obj->sc_frac = st->sc_frac.p;
  obj->cns = st->cns;
  obj->QQ = & st->QQ_p[0];
  obj->AA = & st->AA_p[0];
  obj->QQ_eig = & st->QQ_eig_p[0];
  obj->QQ_imag = & st->QQ_imag_p[0];
  obj->QQ_evec = & st->QQ_evec_p[0];
  obj->QQ_einv = & st->QQ_einv_p[0];
  obj->A = st->A.pp;
  return obj;
}

extern "C" QUBFAST_API void qubx_stim_rates_free(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  delete st;
  delete obj;
}

extern "C" QUBFAST_API void qubx_stim_rates_setup_QQ(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->setup_QQ();
}

extern "C" QUBFAST_API void qubx_stim_rates_setup_AA(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->setup_AA();
}

extern "C" QUBFAST_API void qubx_stim_rates_setup_A(qubx_stim_rates *obj, int sc, double dt)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->setup_A(sc, dt);
}

extern "C" QUBFAST_API void qubx_stim_rates_fix_rate(qubx_stim_rates *obj, int a, int b)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->fixRate(a, b);
}

extern "C" QUBFAST_API void qubx_stim_rates_fix_exp(qubx_stim_rates *obj, int a, int b)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->fixExp(a, b);
}

extern "C" QUBFAST_API void qubx_stim_rates_scale_rate(qubx_stim_rates *obj, int a, int b, int c, int d)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->scaleRate(a, b, c, d);
}

extern "C" QUBFAST_API void qubx_stim_rates_scale_exp(qubx_stim_rates *obj, int a, int b, int c, int d)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->scaleExp(a, b, c, d);
}

extern "C" QUBFAST_API void qubx_stim_rates_balance_loop(qubx_stim_rates *obj, int Nloop, int *loopstates)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->balanceLoop(Nloop, loopstates);
}

extern "C" QUBFAST_API void qubx_stim_rates_balance_loops(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->balanceLoops();
}

extern "C" QUBFAST_API int qubx_stim_rates_setup_constraints(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  return st->setup_constraints();
}

extern "C" QUBFAST_API void qubx_stim_rates_pars_to_model(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->pars_to_model();
}

extern "C" QUBFAST_API void qubx_stim_rates_model_to_pars(qubx_stim_rates *obj)
{
  qubx_stim_rates_storage *st = (qubx_stim_rates_storage *) obj->storage;
  st->model_to_pars();
}

extern "C" QUBFAST_API void qubx_stim_rates_calc_std(qubx_stim_rates *obj, double *InvHessian)
{
  std::vector<double> dPar(obj->cns->Npar);
  qubx_constrained_calc_std(obj->cns, InvHessian, & dPar[0]);
  dK_to_dK012(obj->model->Nstate, obj->cns->pars, &dPar[0], obj->model->K0, obj->model->d_K0, obj->model->d_K1, obj->model->d_K2);
}