<commit_before>/*
 *  DESPOT_Functors.h
 *
 *  Created by <NAME> on 16/08/2012.
 *  Copyright (c) 2012-2013 <NAME>.
 *
 *  This Source Code Form is subject to the terms of the Mozilla Public
 *  License, v. 2.0. If a copy of the MPL was not distributed with this
 *  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 */

#ifndef DESPOT_Functors_h
#define DESPOT_Functors_h

#include <vector>
#include <array>
#include <map>
#include <iostream>
#include <exception>
#include <string>
#include <memory>
#include <Eigen/Dense>

#include "DESPOT.h"

using namespace std;
using namespace Eigen;

//******************************************************************************
#pragma mark Signal Functors
//******************************************************************************
enum class Components {
	One, Two, Three
};

class SignalFunctor {
	public:
		Components m_nC;
		double m_TR;
		ArrayXd m_flip;
		SignalFunctor(const ArrayXd &flip, const double TR, const Components nC);
		virtual ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const = 0;
		virtual size_t size() const { return m_flip.rows(); }
};

class SPGR_Functor : public SignalFunctor {
	public:
		SPGR_Functor(const ArrayXd &flip, const double TR, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
class SPGR_Finite_Functor : public SignalFunctor {
	public:
		double m_Trf, m_TE;
		SPGR_Finite_Functor(const ArrayXd &flip, const double TR, const double Trf, const double TE, const Components nc);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
class SSFP_Functor : public SignalFunctor {
	public:
		ArrayXd m_phases;
		SSFP_Functor(const ArrayXd &flip, const double TR, const ArrayXd &phases, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
		size_t size() const override;
};
class SSFP_Finite_Functor : public SSFP_Functor {
	public:
		double m_Trf, m_TE;
		SSFP_Finite_Functor(const ArrayXd &flip, const double TR, const double Trf, const ArrayXd &phases, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
		
shared_ptr<SignalFunctor> parseSPGR(const Nifti &img, const bool prompt, const Components nC, const bool finite);
shared_ptr<SignalFunctor> parseSSFP(const Nifti &img, const bool prompt, const Components nC, const bool finite);

//******************************************************************************
#pragma mark Optimisation Functor Base Class
//******************************************************************************
// From Nonlinear Tests in Eigen 
template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
class Functor {
	public:
		typedef _Scalar Scalar;
		enum {
			InputsAtCompileTime = NX,
			ValuesAtCompileTime = NY
		};
		typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
		typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
		typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
		
		const long m_inputs, m_values;
		
		Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
		Functor(long inputs, long values) : m_inputs(inputs), m_values(values) {}
		Functor(long values) : m_inputs(InputsAtCompileTime), m_values(values) {}
		
		virtual ~Functor() {};
		
		virtual const long inputs() const { return m_inputs; }
		virtual const long values() const { return m_values; }
		
		virtual int operator()(const Ref<VectorXd> &params, Ref<ArrayXd> diffs) = 0;
		
		virtual const ArrayXd theory(const Ref<VectorXd> &params) = 0;
		virtual const ArrayXd actual() const = 0;
};

//******************************************************************************
#pragma mark Basic DESPOT Functor
//******************************************************************************
class DESPOTFunctor : public Functor<double> {
	public:
		enum class FieldStrength {
			Three, Seven, Unknown
		};
		static const string to_string(const FieldStrength& f) {
			switch (f) {
				case FieldStrength::Three: return "3";
				case FieldStrength::Seven: return "7";
				case FieldStrength::Unknown: return "User";
			}
		}
		enum class Scaling {
			Global, NormToMean
		};
		static const string to_string(const Scaling &p) {
			switch (p) {
				case Scaling::Global: return "global";
				case Scaling::NormToMean: return "normalised to per signal mean";
			}
		}
		enum class OffResMode {
			Map = 0, Single, SingleSymmetric
		};
	
	protected:
		const FieldStrength m_fieldStrength;
		const OffResMode m_offRes;
		const Scaling m_scaling;
		size_t m_nV;
		vector<shared_ptr<SignalFunctor>> m_signals;
		vector<ArrayXd> m_actual, m_theory;
		vector<string> m_names; // Subclasses responsible for initialising this
		const bool m_debug;
	
		ArrayXXd offResBounds() {
			double minTR = numeric_limits<double>::max();
			for (auto &s : m_signals) {
				if (s->m_TR < minTR)
					minTR = s->m_TR;
			}
			switch (m_offRes) {
				case OffResMode::Map: return ArrayXXd(0, 0); break;
				case OffResMode::Single: {
					ArrayXXd b(1, 2);
					b << -0.5 / minTR, 0.5 / minTR;
					return b;
				} break;
				case OffResMode::SingleSymmetric: {
					ArrayXXd b(1, 2);
					b << 0., 0.5 / minTR;
					return b;
				} break;
			}
		}
	
		ArrayXXd PDBounds() {
			ArrayXXd b(nPD(), 2);
			for (size_t i = 0; i < nPD(); i++) {
				b(i, 0) = 1.e4;
				b(i, 1) = 5.e6;
			}
			return b;
		}
		
		virtual const ArrayXXd defaultBounds() = 0;
		
	public:
		double m_f0, m_B1;
		const size_t nOffRes() const {
			switch (m_offRes) {
				case OffResMode::Map: return 0;
				case OffResMode::Single: return 1;
				case OffResMode::SingleSymmetric: return 1;
			}
		}
		
		const size_t nPD() const {
			switch (m_scaling) {
				case Scaling::Global:     return 1;
				case Scaling::NormToMean: return 0;
			}
		}
		
		virtual const size_t nP() const = 0;
		virtual const bool constraint(const VectorXd &params) const = 0;
		const long inputs() const override { return this->nP() + nOffRes() + nPD(); }
		const long values() const override { return m_nV; }
		
		DESPOTFunctor(vector<shared_ptr<SignalFunctor>> &signals_in,
					  const FieldStrength &tesla,
					  const OffResMode &offRes = OffResMode::Single,
					  const Scaling &s = Scaling::NormToMean,
		              const bool &debug = false) :
			m_signals(signals_in),
			m_fieldStrength(tesla), m_offRes(offRes),
			m_scaling(s), m_debug(debug)
		{
			m_actual.reserve(m_signals.size());
			m_theory.reserve(m_signals.size());
			m_nV = 0;
			for (size_t i = 0; i < m_signals.size(); i++) {
				m_actual.emplace_back(m_signals.at(i)->size());
				m_theory.emplace_back(m_signals.at(i)->size());
				m_nV += m_signals.at(i)->size();
			}
		}
		
		const vector<string> &names() { return m_names; }
		shared_ptr<SignalFunctor> &signal(const size_t i) { return m_signals.at(i); }
		ArrayXd &actual(const size_t s) { return m_actual.at(s); }
		const ArrayXd actual() const {
			ArrayXd v(values());
			int index = 0;
			if (m_debug) cout << __PRETTY_FUNCTION__ << endl;
			for (size_t i = 0; i < m_actual.size(); i++) {
				v.segment(index, m_actual.at(i).size()) = m_actual.at(i);
				index += m_actual.at(i).size();
				if (m_debug) cout << v.transpose() << endl;
			}
			return v;
		}
		
		int operator()(const Ref<VectorXd> &params, Ref<ArrayXd> diffs) override {
			eigen_assert(diffs.size() == values());
			ArrayXd t = theory(params);
			ArrayXd s = actual();
			diffs = t - s;
			if (m_debug) {
				cout << endl << __PRETTY_FUNCTION__ << endl;
				cout << "p:      " << params.transpose() << endl;
				cout << "t:      " << t.transpose() << endl;
				cout << "s:      " << s.transpose() << endl;
				cout << "Diffs:  " << diffs.transpose() << endl;
				cout << "Sum:    " << diffs.square().sum() << endl;
			}
			return 0;
		}
};

//******************************************************************************
#pragma mark mcDESPOT Functor
//******************************************************************************
class mcDESPOT : public DESPOTFunctor {
	public:
		static const string to_string(const Components& c) {
			switch (c) {
				case Components::One: return "1";
				case Components::Two: return "2";
				case Components::Three: return "3";
			}
		};
	protected:
		const Components m_components;
		const bool m_finite;
	
	public:
		const size_t nP() const override {
			size_t n = 0;
			switch (m_components) {
				case Components::One: n = 2; break;
				case Components::Two: n = 6; break;
				case Components::Three: n = 9; break;
			}
			if (m_finite)
				n += 1;
			return n;
		}
		
		const ArrayXXd defaultBounds() override {
			ArrayXXd b(inputs(), 2);
			switch (m_fieldStrength) {
				case FieldStrength::Three:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2) << 0.1, 4.0, 0.01, 1.5; break;
						case Components::Two:   b.block(0, 0, 6, 2) << 0.1, 0.25, 0.002, 0.03, 0.7, 4.0, 0.075, 0.145, 0.05, 0.3, 0.0, 0.95; break;
						case Components::Three: b.block(0, 0, 9, 2) << 0.1, 0.25, 0.002, 0.03, 0.7, 2.0, 0.075, 0.145, 3.5, 4.0, 0.8, 1.5, 0.05, 0.3, 0.001, 0.3, 0.0, 0.95; break;
					} break;
				case FieldStrength::Seven:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2) << 0.1, 4.0, 0.01, 2.0; break;
						case Components::Two:   b.block(0, 0, 6, 2) << 0.1, 0.5, 0.001, 0.025, 1.0, 4.0, 0.04, 0.08, 0.01, 0.25, 0.001, 1.0; break;
						case Components::Three: b.block(0, 0, 9, 2) << 0.1, 0.5, 0.001, 0.025, 1.0, 2.5, 0.04, 0.08, 3., 4.5, 0.5, 2.0, 0.01, 0.25, 0.001, 0.4, 0.001, 1.0; break;
					} break;
				case FieldStrength::Unknown:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2).setZero(); break;
						case Components::Two:   b.block(0, 0, 6, 2).setZero(); break;
						case Components::Three: b.block(0, 0, 9, 2).setZero(); break;
					} break;
			}
			if (m_finite) {
				b.block(nP() - 1, 0, 1, 2) << 0., 100.;
			}
			b.block(nP(), 0, nOffRes(), 2) = offResBounds();
			b.block(nP() + nOffRes(), 0, nPD(), 2) = PDBounds();
			return b;
		}
		
		const ArrayXd defaultThresholds() {
			ArrayXd m(inputs());
			switch (m_components) {
				case Components::One: m.head(nP()) << 0.05, 0.05; break;
				case Components::Two: m.head(nP()) << 0.5, 0.5, 0.5, 0.5, 0.5, 0.05; break;
				case Components::Three: m.head(nP()) << 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.05, 0.05; break;
			}<commit_msg>Updated the default 3T bounds (again).<commit_after>/*
 *  DESPOT_Functors.h
 *
 *  Created by <NAME> on 16/08/2012.
 *  Copyright (c) 2012-2013 <NAME>.
 *
 *  This Source Code Form is subject to the terms of the Mozilla Public
 *  License, v. 2.0. If a copy of the MPL was not distributed with this
 *  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 */

#ifndef DESPOT_Functors_h
#define DESPOT_Functors_h

#include <vector>
#include <array>
#include <map>
#include <iostream>
#include <exception>
#include <string>
#include <memory>
#include <Eigen/Dense>

#include "DESPOT.h"

using namespace std;
using namespace Eigen;

//******************************************************************************
#pragma mark Signal Functors
//******************************************************************************
enum class Components {
	One, Two, Three
};

class SignalFunctor {
	public:
		Components m_nC;
		double m_TR;
		ArrayXd m_flip;
		SignalFunctor(const ArrayXd &flip, const double TR, const Components nC);
		virtual ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const = 0;
		virtual size_t size() const { return m_flip.rows(); }
};

class SPGR_Functor : public SignalFunctor {
	public:
		SPGR_Functor(const ArrayXd &flip, const double TR, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
class SPGR_Finite_Functor : public SignalFunctor {
	public:
		double m_Trf, m_TE;
		SPGR_Finite_Functor(const ArrayXd &flip, const double TR, const double Trf, const double TE, const Components nc);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
class SSFP_Functor : public SignalFunctor {
	public:
		ArrayXd m_phases;
		SSFP_Functor(const ArrayXd &flip, const double TR, const ArrayXd &phases, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
		size_t size() const override;
};
class SSFP_Finite_Functor : public SSFP_Functor {
	public:
		double m_Trf, m_TE;
		SSFP_Finite_Functor(const ArrayXd &flip, const double TR, const double Trf, const ArrayXd &phases, const Components nC);
		ArrayXd signal(const VectorXd &p, const double B1 = 1., const double f0 = 0.) const override;
};
		
shared_ptr<SignalFunctor> parseSPGR(const Nifti &img, const bool prompt, const Components nC, const bool finite);
shared_ptr<SignalFunctor> parseSSFP(const Nifti &img, const bool prompt, const Components nC, const bool finite);

//******************************************************************************
#pragma mark Optimisation Functor Base Class
//******************************************************************************
// From Nonlinear Tests in Eigen 
template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
class Functor {
	public:
		typedef _Scalar Scalar;
		enum {
			InputsAtCompileTime = NX,
			ValuesAtCompileTime = NY
		};
		typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
		typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
		typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
		
		const long m_inputs, m_values;
		
		Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
		Functor(long inputs, long values) : m_inputs(inputs), m_values(values) {}
		Functor(long values) : m_inputs(InputsAtCompileTime), m_values(values) {}
		
		virtual ~Functor() {};
		
		virtual const long inputs() const { return m_inputs; }
		virtual const long values() const { return m_values; }
		
		virtual int operator()(const Ref<VectorXd> &params, Ref<ArrayXd> diffs) = 0;
		
		virtual const ArrayXd theory(const Ref<VectorXd> &params) = 0;
		virtual const ArrayXd actual() const = 0;
};

//******************************************************************************
#pragma mark Basic DESPOT Functor
//******************************************************************************
class DESPOTFunctor : public Functor<double> {
	public:
		enum class FieldStrength {
			Three, Seven, Unknown
		};
		static const string to_string(const FieldStrength& f) {
			switch (f) {
				case FieldStrength::Three: return "3";
				case FieldStrength::Seven: return "7";
				case FieldStrength::Unknown: return "User";
			}
		}
		enum class Scaling {
			Global, NormToMean
		};
		static const string to_string(const Scaling &p) {
			switch (p) {
				case Scaling::Global: return "global";
				case Scaling::NormToMean: return "normalised to per signal mean";
			}
		}
		enum class OffResMode {
			Map = 0, Single, SingleSymmetric
		};
	
	protected:
		const FieldStrength m_fieldStrength;
		const OffResMode m_offRes;
		const Scaling m_scaling;
		size_t m_nV;
		vector<shared_ptr<SignalFunctor>> m_signals;
		vector<ArrayXd> m_actual, m_theory;
		vector<string> m_names; // Subclasses responsible for initialising this
		const bool m_debug;
	
		ArrayXXd offResBounds() {
			double minTR = numeric_limits<double>::max();
			for (auto &s : m_signals) {
				if (s->m_TR < minTR)
					minTR = s->m_TR;
			}
			switch (m_offRes) {
				case OffResMode::Map: return ArrayXXd(0, 0); break;
				case OffResMode::Single: {
					ArrayXXd b(1, 2);
					b << -0.5 / minTR, 0.5 / minTR;
					return b;
				} break;
				case OffResMode::SingleSymmetric: {
					ArrayXXd b(1, 2);
					b << 0., 0.5 / minTR;
					return b;
				} break;
			}
		}
	
		ArrayXXd PDBounds() {
			ArrayXXd b(nPD(), 2);
			for (size_t i = 0; i < nPD(); i++) {
				b(i, 0) = 1.e4;
				b(i, 1) = 5.e6;
			}
			return b;
		}
		
		virtual const ArrayXXd defaultBounds() = 0;
		
	public:
		double m_f0, m_B1;
		const size_t nOffRes() const {
			switch (m_offRes) {
				case OffResMode::Map: return 0;
				case OffResMode::Single: return 1;
				case OffResMode::SingleSymmetric: return 1;
			}
		}
		
		const size_t nPD() const {
			switch (m_scaling) {
				case Scaling::Global:     return 1;
				case Scaling::NormToMean: return 0;
			}
		}
		
		virtual const size_t nP() const = 0;
		virtual const bool constraint(const VectorXd &params) const = 0;
		const long inputs() const override { return this->nP() + nOffRes() + nPD(); }
		const long values() const override { return m_nV; }
		
		DESPOTFunctor(vector<shared_ptr<SignalFunctor>> &signals_in,
					  const FieldStrength &tesla,
					  const OffResMode &offRes = OffResMode::Single,
					  const Scaling &s = Scaling::NormToMean,
		              const bool &debug = false) :
			m_signals(signals_in),
			m_fieldStrength(tesla), m_offRes(offRes),
			m_scaling(s), m_debug(debug)
		{
			m_actual.reserve(m_signals.size());
			m_theory.reserve(m_signals.size());
			m_nV = 0;
			for (size_t i = 0; i < m_signals.size(); i++) {
				m_actual.emplace_back(m_signals.at(i)->size());
				m_theory.emplace_back(m_signals.at(i)->size());
				m_nV += m_signals.at(i)->size();
			}
		}
		
		const vector<string> &names() { return m_names; }
		shared_ptr<SignalFunctor> &signal(const size_t i) { return m_signals.at(i); }
		ArrayXd &actual(const size_t s) { return m_actual.at(s); }
		const ArrayXd actual() const {
			ArrayXd v(values());
			int index = 0;
			if (m_debug) cout << __PRETTY_FUNCTION__ << endl;
			for (size_t i = 0; i < m_actual.size(); i++) {
				v.segment(index, m_actual.at(i).size()) = m_actual.at(i);
				index += m_actual.at(i).size();
				if (m_debug) cout << v.transpose() << endl;
			}
			return v;
		}
		
		int operator()(const Ref<VectorXd> &params, Ref<ArrayXd> diffs) override {
			eigen_assert(diffs.size() == values());
			ArrayXd t = theory(params);
			ArrayXd s = actual();
			diffs = t - s;
			if (m_debug) {
				cout << endl << __PRETTY_FUNCTION__ << endl;
				cout << "p:      " << params.transpose() << endl;
				cout << "t:      " << t.transpose() << endl;
				cout << "s:      " << s.transpose() << endl;
				cout << "Diffs:  " << diffs.transpose() << endl;
				cout << "Sum:    " << diffs.square().sum() << endl;
			}
			return 0;
		}
};

//******************************************************************************
#pragma mark mcDESPOT Functor
//******************************************************************************
class mcDESPOT : public DESPOTFunctor {
	public:
		static const string to_string(const Components& c) {
			switch (c) {
				case Components::One: return "1";
				case Components::Two: return "2";
				case Components::Three: return "3";
			}
		};
	protected:
		const Components m_components;
		const bool m_finite;
	
	public:
		const size_t nP() const override {
			size_t n = 0;
			switch (m_components) {
				case Components::One: n = 2; break;
				case Components::Two: n = 6; break;
				case Components::Three: n = 9; break;
			}
			if (m_finite)
				n += 1;
			return n;
		}
		
		const ArrayXXd defaultBounds() override {
			ArrayXXd b(inputs(), 2);
			switch (m_fieldStrength) {
				case FieldStrength::Three:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2) << 0.1, 4.0, 0.01, 1.5; break;
						case Components::Two:   b.block(0, 0, 6, 2) << 0.200, 0.350, 0.005, 0.015, 0.700, 2.000, 0.050, 0.120, 0.050, 0.200, 0.0, 0.5; break;
						case Components::Three: b.block(0, 0, 9, 2) << 0.200, 0.350, 0.005, 0.015, 0.700, 2.000, 0.050, 0.120, 3.500, 7.000, 3.000, 7.000, 0.050, 0.200, 0, 0.5, 0, 1.0; break;
					} break;
				case FieldStrength::Seven:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2) << 0.1, 4.0, 0.01, 2.0; break;
						case Components::Two:   b.block(0, 0, 6, 2) << 0.1, 0.5, 0.001, 0.025, 1.0, 4.0, 0.04, 0.08, 0.01, 0.25, 0.001, 1.0; break;
						case Components::Three: b.block(0, 0, 9, 2) << 0.1, 0.5, 0.001, 0.025, 1.0, 2.5, 0.04, 0.08, 3., 4.5, 0.5, 2.0, 0.01, 0.25, 0.001, 0.4, 0.001, 1.0; break;
					} break;
				case FieldStrength::Unknown:
					switch (m_components) {
						case Components::One:   b.block(0, 0, 2, 2).setZero(); break;
						case Components::Two:   b.block(0, 0, 6, 2).setZero(); break;
						case Components::Three: b.block(0, 0, 9, 2).setZero(); break;
					} break;
			}
			if (m_finite) {
				b.block(nP() - 1, 0, 1, 2) << 0., 100.;
			}
			b.block(nP(), 0, nOffRes(), 2) = offResBounds();
			b.block(nP() + nOffRes(), 0, nPD(), 2) = PDBounds();
			return b;
		}
		
		const ArrayXd defaultThresholds() {
			ArrayXd m(inputs());
			switch (m_components) {
				case Components::One: m.head(nP()) << 0.05, 0.05; break;
				case Components::Two: m.head(nP()) << 0.5, 0.5, 0.5, 0.5, 0.5, 0.05; break;
				case Components::Three: m.head(nP()) << 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.05, 0.05; break;
			}
{'fraction_non_alphanumeric': 0.19423231987857514, 'fraction_numerical': 0.03552625337950007, 'mean_word_length': 3.618619934282585, 'pattern_counts': {'":': 0, '<': 209, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5485717', 'n_tokens_mistral': 9110, 'n_tokens_neox': 7435, 'n_words': 2055}