/*
 Copyright (C) 2022 Leiden University Medical Center

 This program 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.

 This program 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
 along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <APqr7.h>
#include <math.h>
#include <vector>

/*
 *********
 * APqr7 *
 *********

This software provides an Action Potential Cure (APqr) to correct divergent
membrane potentials in excitable biological systems. The first X action
potentials are logged when the software starts, after which AP correction
starts from the (X+1)-st AP onwards. This correction occurs with the use of
dynamic pacth-clamp currrent injection.

IN:
	*) Cm				Capacitance of the cell
	*) V_cutoff			Threshold potential for the detection of the beginning
						of an AP
	*) Slope_tresh		Slope threshold that defines the beginning of the
						AP (mV/ms)
	*) BCL_cutoff		Threshold value for the end of an AP, given as a
						percentage of the total APD
	*) lognum			Number of APs that need to be logged as a reference
	*) Rm				Initial resistance
	*) Rm_corr_up		Factor to increase Rm with when necessary
	*) Rm_corr_down		Factor to decrease Rm with when necessary
	*) noise_tresh		The noise level that is allowed around the ideal value
						before correcting
OUT:
	*) Vout 			voltage that is used to inject the calculated amount
						of current into the excitable system
*/

/*
createRTXIPlugin
----------------
Creation of a new RTXI Plugin

IN:
	*) None
OUT:
	*) RTXIPlugin
*/
extern "C" Plugin::Object *createRTXIPlugin(void)
{
	return new gAPqr7();
}

/*
vars[]
----
This is not a function, but rather the construction of a list. This list contains
all the variables that are visible and/or modifiable in the GUI of the software
module. There is the choice between INPUT, OUTPUT, PARAMETER and STATE.

INPUT: connected to the input port
OUTPUT: connected to the output port
PARAMETER: modifiable variable in the code
STATE: non-modifiable variable in the code
*/
static DefaultGUIModel::variable_t vars[] = {
	{ "Vm (mV)", "Membrane potential (mV)", DefaultGUIModel::INPUT, },
	{ "Iout (pA)", "Output current (pA)", DefaultGUIModel::OUTPUT, },
	{ "iAP", "ideal AP", DefaultGUIModel::OUTPUT, },
	{ "Cm (pF)", "pF", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, },
	{ "V_cutoff (mV)", "Threshold potential for the detection of the beginning of an AP, together with Slope_thresh",
	DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
	{ "Slope_thresh (mV/ms)", "Slope threshold that defines the beginning of the AP (mV/ms)",
	DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
	{ "BCL_cutoff (pct)", "Threshold value for the end of an AP, given as a percentage of the total APD",
	DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
	{ "noise_tresh (mV)", "The noise level that is allowed before correcting", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, }, 
	{ "Rm (MOhm)", "MOhm", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, },
	{ "lognum", "Number of APs that need to be logged as a reference", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, },
	{ "Rm_corr_up", "To increase Rm when necessary", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, },
	{ "Rm_corr_down", "To decrease Rm when necessary", DefaultGUIModel::PARAMETER
	| DefaultGUIModel::DOUBLE, }, 
	{ "Correction (0 or 1)", "Switch Rm correction off (0) or on (1)",
	DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
	{ "Vm2 (mV)", "Membrane potential (mV)", DefaultGUIModel::STATE, },
	{ "Iout2 (pA)", "Output Current (pA)", DefaultGUIModel::STATE, },
	{ "Period (ms)", "Period (ms)", DefaultGUIModel::STATE, }, 
	{ "Time (ms)", "Time (ms)", DefaultGUIModel::STATE, },
	{ "APs2", "APs", DefaultGUIModel::STATE, },
	{ "log_ideal_on2", "log_ideal_on", DefaultGUIModel::STATE, },
	{ "BCL2", "BCL", DefaultGUIModel::STATE, },
	{ "enter2", "enter", DefaultGUIModel::STATE, },
	{ "Rm2 (MOhm)", "MOhm", DefaultGUIModel::STATE, },
	{ "act2", "0 or 1", DefaultGUIModel::STATE, },
	{ "count", "number", DefaultGUIModel::STATE, },
	{ "count2", "number", DefaultGUIModel::STATE, },
	{ "modulo_state", "number", DefaultGUIModel::STATE, },
};

/*
num_vars
--------

*/
static size_t num_vars = sizeof(vars) / sizeof(DefaultGUIModel::variable_t);

/*
gAPqr7
------

IN:
	*)
OUT:
	*)
*/
gAPqr7::gAPqr7(void) : DefaultGUIModel("APqr7", ::vars, ::num_vars)
{
	setWhatsThis(
		"<p><b>APqr:</b><br>APqr7 </p>");
	DefaultGUIModel::createGUI(vars, num_vars);
	initParameters();
	update(INIT);
	refresh();
	resizeMe();
}

gAPqr7::~gAPqr7(void) {}

/*
cleanup
-------

IN:
	*) None
OUT:
	*) None
*/
void gAPqr7::cleanup()
{
	for(i=0;i<10000;i++){
		Vm_log[i]=0;
		Vm_diff_log[i]=0;
		ideal_AP[i]=0;
	}
}

/*
execute
-------

IN:
	*) None
OUT:
	*) Vout				voltage that is used to inject the calculated amount
						of current into the excitable system
*/
void gAPqr7::execute(void)
{
	systime = count * period; // time in milli-seconds
	Vm = input(0) * 1e2; // convert 10V to mV. Divided by 10 because the amplifier produces 10-fold amplified voltages. Multiplied by 1000 to convert V to mV.

	Vm_log[count % (int)modulo] = Vm;

	if(count>(int)(1/period)-1 && (Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) >= slope_thresh && APs<lognum && enter == 0 && Vm > V_cutoff)
	{
		BCL = (APs==-1? 0: (BCL*APs + count2)/(APs+1));
		log_ideal_on = 1;
		count2 = 0;
		enter = 1;
		APs++;
	}

	if((Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) < 0 && enter == 1)
	{
		enter = 0;
	}

	if(APs<lognum && log_ideal_on == 1)
	{
		ideal_AP[count2] = (ideal_AP[count2]*APs + Vm)/(APs+1);
		count2++;
	}

	if (act == 0 && (Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) >= slope_thresh && APs >= lognum && Vm > V_cutoff)
	{
		count = 0;		
		act = 1;
	}

	
	if (act == 1)
	{
		Iout = Cm * (1/Rm) * (Vm - ideal_AP[count]);
		output(0) = -Iout * 2.5e-3; // This is equal to Vout
		Vm_diff_log[count] = Vm - ideal_AP[count];

		iAP = ideal_AP[count];
		output(1) = iAP;
	}

	if(corr == 1 && act == 1 && count > 1 && abs(Vm_diff_log[count])>noise_tresh)
	{
		if((Vm_diff_log[count-1] / Vm_diff_log[count]) < 0) //they have the opposite sign, so an overshoot occured
		{
			Rm = Rm * Rm_corr_up;
		}
		if(abs(Vm_diff_log[count-1]) < abs(Vm_diff_log[count]) && (Vm_diff_log[count-1] / Vm_diff_log[count]) > 0) //the difference is getting larger, so increase current
		{
			Rm = Rm / Rm_corr_down;
		}
	}


	if (count > BCL_cutoff*BCL)
	{
		act = 0;
		output(0) = 0;
	}

	count++;

	count_r = (double)count/1000.0;
	count2_r = (double)count2/1000.0;

}

/*
Update
------
ABC

IN:
	*) flag				Indicating the state of the update:
						INIT, MODIFY, PERIOD, PAUSE, UNPAUSE
OUT:
	*) None
*/
void gAPqr7::update(DefaultGUIModel::update_flags_t flag)
{
	switch (flag)
	{
	case INIT:
		setParameter("Cm (pF)", Cm);
		setParameter("V_cutoff (mV)", V_cutoff);
		setParameter("Rm (MOhm)", Rm);
		setParameter("Rm_corr_up", Rm_corr_up);
		setParameter("Rm_corr_down", Rm_corr_down);
		setParameter("noise_tresh (mV)", noise_tresh);
		setParameter("lognum", lognum);
		setParameter("BCL_cutoff (pct)", BCL_cutoff);
		setParameter("Slope_thresh (mV/ms)", slope_thresh);
		setParameter("Correction (0 or 1)", corr);
		setState("Vm2 (mV)", Vm);
		setState("Iout2 (pA)", Iout);
		setState("Time (ms)", systime);
		setState("Period (ms)", period);
		setState("APs2", APs);
		setState("log_ideal_on2", log_ideal_on);
		setState("BCL2", BCL);
		setState("enter2", enter);
		setState("Rm2 (MOhm)", Rm);
		setState("act2", act);
		setState("count", count_r);
		setState("count2", count2_r);
		setState("modulo_state", modulo);
		break;
	case MODIFY:
		Cm = getParameter("Cm (pF)").toDouble();
		Rm = getParameter("Rm (MOhm)").toDouble();
		lognum = getParameter("lognum").toDouble();
		V_cutoff = getParameter("V_cutoff (mV)").toDouble();
		BCL_cutoff = getParameter("BCL_cutoff (pct)").toDouble();
		noise_tresh = getParameter("noise_tresh (mV)").toDouble();
		Rm_corr_up = getParameter("Rm_corr_up").toDouble();
		Rm_corr_down = getParameter("Rm_corr_down").toDouble();
		slope_thresh = getParameter("Slope_thresh (mV/ms)").toDouble();
		corr = getParameter("Correction (0 or 1)").toDouble();		
		systime = 0;
		count = 0;
		APs = -1;
		BCL = 0;
		log_ideal_on = 0;
		enter = 0;
		count2 = 0;
		cleanup();
		break;
	case PERIOD:
		period = RT::System::getInstance()->getPeriod() * 1e-6; // time in milli-seconds
		modulo = (1.0/(RT::System::getInstance()->getPeriod() * 1e-6)) * 1000.0;
		break;
	case PAUSE:
		output(0) = 0.0;
		Iout = 0;
		act = 0;
		systime = 0;
		break;
	case UNPAUSE:
		break;
	default:
		break;
	}
}

/*
initParameters
--------------
This function sets all values to their defaults when no external parameters are provided
through the GUI interface.

IN:
	*) None
OUT:
	*) None
*/
void gAPqr7::initParameters()
{
	Vm = -80; 			// mV
	Cm = 150; 			// pF
	Rm = 150; 			// MOhm
	slope_thresh = 5.0; // mV
	corr = 1;
	Iout = 0;			// pA
	output(0) = -Iout * 0.5e-3;
	period = RT::System::getInstance()->getPeriod() * 1e-6; // ms
	systime = 0;
	count = 0;
	act = 0;
	iAP=0;
	Rm_corr_up=8;
	Rm_corr_down=2;
	noise_tresh = 0.5; 	// mV
	BCL = 0;			// ms
	count2 = 0;
	APs = -1;
	V_cutoff = -40;		// mV
	BCL_cutoff = 0.98;
	enter = 0;
	log_ideal_on = 0;
	lognum = 3;
	count_r = 0;
	count2_r = 0;
	modulo = (1.0/(RT::System::getInstance()->getPeriod() * 1e-6)) * 1000.0;
}