$$\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 QUBX_MODEL_H
#define QUBX_MODEL_H

#include "qubfast.h"

#ifdef __cplusplus
extern "C" {
#endif

// Model has constant rates, amps etc.
#pragma pack(push, 1) 
typedef struct {
  // maintained by model:
  int version;
  void *storage;
  int Nstate, Nclass;
  int *clazz; // clazz[i] = class_index of state i
  // you set these:
  int usePeq, IeqFv;
  int Nchannel;
  int iVoltage; // 1-based index of voltage signal (since 0 is the current)
#ifdef EMSCRIPTEN
  int r1; // emscripten wants doubles aligned at 64
#endif
  double Gleak, Vleak, Vrev;
  // and fill these:
  double amp0, var0; // per-class; honored only if IeqFv
  double *amp, *var; // per class; units of IeqFv ? pS : pA;  IeqFv: excess conductance;  else: total single-channel current
  double *P0; // per state
  double **K0, **K1;
  int **L, **V;
  
  // set nonzero to interrupt
  int stop_flag;
  int progress; // 0 to 100

#ifdef EMSCRIPTEN
  int r2;
#endif
  // see calc_std funcs:
  double d_Nchannel;
  double d_amp0, d_var0;
  double *d_amp, *d_var;
  double **d_K0, **d_K1;

  // end version 1
  
  double **K2, **d_K2;
  int **P;
  // end version 2
} qubx_model;
#pragma pack(pop) 

QUBFAST_API qubx_model* qubx_model_create(int Nstate, int Nclass, int *class_of_state,
					      callbk_reportfun report_cb, void *report_data);
QUBFAST_API void qubx_model_free(qubx_model *m);


typedef struct {
  int version;
  void *storage;
  qubx_model *single_model;
  qubx_model *multi_model;
  int nmutistate, nmutipath, nmuticlass, nclass_order, npath;
  int *class_order, *cmutistate, *path, *mutistate, *mutipath;
  // end version 1
} qubx_multi_model;

QUBFAST_API void qubx_mutimakv_size(int nchannel, int nstate, int npath, int *path,
				    int *nmutistate, int *nmutipath);
// returns upper bound on array sizes needed for mutistate and mutipath, via args

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]


QUBFAST_API qubx_multi_model* qubx_multi_model_create(qubx_model *single_model);
QUBFAST_API void qubx_multi_model_free(qubx_multi_model *mm);
QUBFAST_API void qubx_multi_model_update(qubx_multi_model *mm); // single_model -> multi_model amps and rates



typedef struct {
  int version;
  void *storage;
  int Npar, NparR, NfPar, Ncns, NcnsR;
  double **Ain, *Bin;
  double *pars, *fPars;
  double **Acns, **Ascal;
  double *Bcns, *Bscal;
  // end version 1
} qubx_constrained;
  
QUBFAST_API qubx_constrained* qubx_constrained_create(int Npar);
QUBFAST_API void qubx_constrained_free(qubx_constrained *obj);
QUBFAST_API int qubx_constrained_setup(qubx_constrained *obj);
QUBFAST_API void qubx_constrained_free_to_pars(qubx_constrained *obj);
QUBFAST_API void qubx_constrained_pars_to_free(qubx_constrained *obj);
QUBFAST_API void qubx_constrained_calc_std(qubx_constrained *obj, double *InvHessian, double *dPar);


typedef struct {
  int version;
  void *storage;
  qubx_model *model;
  int Nstim, Nsig;
  double *stimclasses;
  double *sc_frac;
  qubx_constrained *cns;
  double **sc_amp, **sc_var, **sc_std;
  // end version 1
} qubx_stim_amps;
// define param order: amp[0], var[0], amp[1], var[1], ..., amp0, var0

QUBFAST_API qubx_stim_amps*  qubx_stim_amps_create(qubx_model *m, int Nstim_, int Nsig_,
						  double *stimclasses_, double *sc_frac_);
QUBFAST_API void qubx_stim_amps_free(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_setup_scs(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_pars_to_model(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_model_to_pars(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_fix_amp0(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_fix_var0(qubx_stim_amps *obj);
QUBFAST_API void qubx_stim_amps_fix_amp(qubx_stim_amps *obj, int cls);
QUBFAST_API void qubx_stim_amps_fix_var(qubx_stim_amps *obj, int cls);
QUBFAST_API void qubx_stim_amps_calc_std(qubx_stim_amps *obj, double *InvHessian);

typedef int (*eigen_func)(int Nstate, double *a, double *wr, double *wi, double *v, double *v_inv);

typedef struct {
  int version;
  void *storage;
  qubx_model *model;
  eigen_func custom_eigen; // you can do better than the built-in (e.g. wrap scipy or linpack)
  int Nstim, Nsig;
  double *stimclasses;
  double *sc_frac;
    
  qubx_constrained *cns;
   
  double ***QQ, ***AA, **QQ_eig, **QQ_imag, ***QQ_evec, ***QQ_einv;
  double **A; // of requested sc, t
  // end version 1
} qubx_stim_rates;
// param order from qubx_constraints: K01toK
  
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);
QUBFAST_API void qubx_stim_rates_free(qubx_stim_rates *obj);
QUBFAST_API void qubx_stim_rates_setup_QQ(qubx_stim_rates *obj); // from rates
QUBFAST_API void qubx_stim_rates_setup_AA(qubx_stim_rates *obj);  // and QQ_*; from QQ
QUBFAST_API void qubx_stim_rates_setup_A(qubx_stim_rates *obj, int sc, double dt);
QUBFAST_API void qubx_stim_rates_fix_rate(qubx_stim_rates *obj, int a, int b);
QUBFAST_API void qubx_stim_rates_fix_exp(qubx_stim_rates *obj, int a, int b);
QUBFAST_API void qubx_stim_rates_scale_rate(qubx_stim_rates *obj, int a, int b, int c, int d);
QUBFAST_API void qubx_stim_rates_scale_exp(qubx_stim_rates *obj, int a, int b, int c, int d);
QUBFAST_API void qubx_stim_rates_balance_loop(qubx_stim_rates *obj, int Nloop, int *loopstates);
QUBFAST_API void qubx_stim_rates_balance_loops(qubx_stim_rates *obj);
QUBFAST_API int qubx_stim_rates_setup_constraints(qubx_stim_rates *obj);
QUBFAST_API void qubx_stim_rates_pars_to_model(qubx_stim_rates *obj);
QUBFAST_API void qubx_stim_rates_model_to_pars(qubx_stim_rates *obj);
QUBFAST_API void qubx_stim_rates_calc_std(qubx_stim_rates *obj, double *InvHessian);


#ifdef __cplusplus
}
#endif

#endif