Fast Automatic Differentiation in C++

//////////////////////////////////////////////////////////////////////////
/// The Reed Library 
/// Copyright (c) 2006 Yevgeniy Shatokhin 

/// File: asection.cpp 
/// Contains implementation of the active section class.

/// 
/// Last update: August 29, 2006

#include "stdafx.h"

#include <reed/adouble.h>
#include <reed/asection.h>

#pragma warning(push)	// store current state for each warning
#pragma warning(disable: 4311 4312)	// ptr-to-int conversion (and vice versa)

namespace reed
{
	//////////////////////////////////////////////////////////////////////////
	// Implementation: CActiveSection
	//////////////////////////////////////////////////////////////////////////
	
	//////////////////////////////////////////////////////////////////////////
	// Static initialization.
	CActiveSection* CActiveSection::activeAS = NULL;

	//////////////////////////////////////////////////////////////////////////
	
	// Constructor.
	CActiveSection::CActiveSection()
		: mode(AS_MODE_PLAIN),
		bCacheEnabled(FALSE)
	{}

	// Copy constructor.
	CActiveSection::CActiveSection(const CActiveSection& rhs)
		: mode(rhs.mode),
		bCacheEnabled(rhs.bCacheEnabled),
		alloc(rhs.alloc)
	{}

	// Destructor.
	CActiveSection::~CActiveSection()
	{}

	// Begin AS.
	void
	CActiveSection::begin()
	{
		// we should be outside any AS
		assert(activeAS == NULL);
		activeAS = this;

		// AS mode is AS_MODE_PLAIN by default. Other modes can be set in end()
		// or within the AS.

		if (mode == AS_MODE_PLAIN)
		{
			// Clear the previously recorded CG  if caching is not used.
			// Any pointers and references to the nodes of this CG become invalid.
			alloc.clear();
		}
		else if (mode == AS_MODE_CACHING1)
		{
			// A node for the current operation is at the head of the CG.
			pCurrent = alloc.getHead();

			// Can be NULL only if there are no nodes in the CG and that could 
			// only result from the user's error. 
			// Anyway, what purpose is there in the empty active section?!
			assert(pCurrent != NULL);	
		}
	}

	// End AS.
	void
	CActiveSection::end()
	{
		// set AS mode for the possible next execution of the same AS
		mode = (bCacheEnabled) ? AS_MODE_CACHING1 : AS_MODE_PLAIN;
		activeAS = NULL;
	}

	int 
	CActiveSection::update(TADVector& arg, TDoubleVector& val)
	{
		assert(arg.size() == val.size());

		TDoubleVector stored(arg.size());	// here the stored previous values reside

		// update specified independent variables, store their previous values
		TADVector::iterator iarg = arg.begin();
		TDoubleVector::iterator ival = val.begin(), 
			istored = stored.begin();

		for (; iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			*istored = iarg->pNode->value.d;
			iarg->pNode->value.d = *ival;

			++ival; ++istored;
		}

		CNodeBase* pHead = alloc.getHead();
		CNodeBase* pStop = alloc.getFree();

		// walk the CG
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			assert(p != NULL);
			if (!p->evaluate()) 
			{
				// The values specified require a CG of different structure.
				// Restore previous values of the independent variables and return.
				istored = stored.begin();
				for (iarg = arg.begin(); iarg != arg.end(); ++iarg)
				{
					iarg->pNode->value.d = *istored;
					++istored;
				}
				return FALSE;
			}
		}

		// OK, update 'value' fields in the nodes
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			p->value.d = p->accum.d;	// boolean values are copied this way too
		}

		return TRUE;
	}

	void
	CActiveSection::getAccumulators(TADVector& arg, TDoubleVector& val)
	{
		val.clear();

		TADVector::iterator iarg = arg.begin();

		for (; iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			val.push_back(iarg->pNode->accum.d);
		}
	}

	void
	CActiveSection::setAccumulators(TADVector& arg, TDoubleVector& val)
	{
		assert(arg.size() == val.size());

		TADVector::iterator iarg = arg.begin();
		TDoubleVector::iterator ival = val.begin();

		for (; iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			iarg->pNode->accum.d = *ival;
			++ival;
		}
	}

	// forward mode
	int
	CActiveSection::forward(TADVector& arg, TDoubleVector& s)
	{
		assert(arg.size() == s.size());
        
		// set accumulators first
		CNodeBase* pHead = alloc.getHead();
		CNodeBase* pStop = alloc.getFree();

		// clear previous accumulated values
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			assert(p != NULL);
			p->accum.d = 0.0;
		}

		TADVector::iterator iarg = arg.begin();
		TDoubleVector::iterator is = s.begin();

        for (; iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			iarg->pNode->accum.d = *is;
			++is;
		}

		int status = TRUE;

        // walk the CG
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			assert(p != NULL);
			int st = p->forwardStep();
			/*
			if (st != TRUE)
						{
							if ((st == FALSE) && (status == TRUE))
							{
								status = FALSE;
							}
							else
							{
								status = FALSE;
							}
						}*/
			status &= st;
		}

		return status;
	}

	int
	CActiveSection::forward(CADouble arg, CADouble func, double& out)
	{
		TADVector av(1);
		av[0] = arg;

		TDoubleVector dv(1);
		dv[0] = 1.0;

		int st = forward(av, dv);
		
		assert(func.pNode);
		out = func.pNode->accum.d;

		return st;
	}

	int
	CActiveSection::forward(TADVector& arg, TADVector& func, TDoubleVector& s, TDoubleVector& out)
	{
		assert(arg.size() == s.size());
		
		int st = forward(arg, s);

		out.clear();
		out.reserve(func.size());
		for (TADVector::iterator ifunc = func.begin(); ifunc != func.end(); ++ifunc)
		{
			assert(ifunc->pNode != NULL);
			out.push_back(ifunc->pNode->accum.d);
		}
		return st;
	}

	// reverse mode
	int 
	CActiveSection::reverse(TADVector& func, TDoubleVector& w)
	{
		assert(func.size() == w.size());

		// set accumulators first
		CNodeBase* pHead = alloc.getHead();
		CNodeBase* pStop = alloc.getFree();
		CNodeBase *p, *next, *prev;

		assert(pHead != NULL);	// at least one node must exist.

		// clear previous accumulated values and prepare 'next' pointers  
		// for the reverse pass (store ptrs to the previous elements in them - XOR'ed)
		prev = NULL;
		next = pHead->next;
		for (p = pHead; next != pStop; p = next)
		{
			assert(p != NULL);
			next = p->next;
			p->next = (CNodeBase*)((unsigned long)next ^ (unsigned long)prev);
			prev = p;
			p->accum.d = 0.0;
		}

		TADVector::iterator ifunc = func.begin();
		TDoubleVector::iterator iw = w.begin();

		for (; ifunc != func.end(); ++ifunc)
		{
			assert(ifunc->pNode != NULL);
			ifunc->pNode->accum.d = *iw;
			++iw;
		}

		int status = TRUE;

		// walk the CG backwards
		// (prev points to the last CG node now)
		while (prev != NULL)
		{
			p = prev;

			// process the current node
			int st = p->reverseStep();
			/*
			if (st != TRUE)
						{
							if ((st == FALSE) && (status == TRUE))
							{
								status = FALSE;
							}
							else
							{
								status = FALSE;
							}
						}*/
			status &= st;

			// restore the 'next' ptr in this node and go to the previous node
			prev = (CNodeBase*)((unsigned long)(p->next) ^ (unsigned long)next);
			p->next = next;
			next = p;
		}
		assert(prev == NULL);	// we must have achieved the head now, prev is 0 for the head.

		return status;
	}

	int 
	CActiveSection::reverse(CADouble arg, CADouble func, double& out)
	{
		TADVector av(1);
		av[0] = func;

		TDoubleVector dv(1);
		dv[0] = 1.0;

		int st = reverse(av, dv);

		assert(arg.pNode);
		out = arg.pNode->accum.d;

		return st;
	}

	int 
	CActiveSection::reverse(TADVector& arg, TADVector& func, TDoubleVector& w, TDoubleVector& out)
	{
		assert(func.size() == w.size());
		
		int st = reverse(func, w);

		out.clear();
		out.reserve(arg.size());
		for (TADVector::iterator iarg = arg.begin(); iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			out.push_back(iarg->pNode->accum.d);
		}
		return st;
	}

	// ree - lower error estimate
	int 
	CActiveSection::ree_lower(TADVector& arg, TDoubleVector& err)
	{
		assert(arg.size() == err.size());

		// set accumulators first
		CNodeBase* pHead = alloc.getHead();
		CNodeBase* pStop = alloc.getFree();

		// set default accumulated values
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			assert(p != NULL);
			// x is assumed to equal x0 <=> x lies in [x0 - 0.5*derr(x0), x0 + 0.5*derr(x0)]
			// It's not difficult to check this even in a debug watch window.
			p->accum.d = 0.5 * derr(p->value.d);
		}

		TADVector::iterator iarg = arg.begin();
		TDoubleVector::iterator ierr = err.begin();

		for (; iarg != arg.end(); ++iarg)
		{
			assert(iarg->pNode != NULL);
			iarg->pNode->accum.d = *ierr;
			++ierr;
		}

		int status = TRUE;

		// walk the CG
		for (CNodeBase* p = pHead; p != pStop; p = p->next)
		{
			assert(p != NULL);
			int st = p->reeStep();
			/*
			if (st != TRUE)
						{
							if ((st == FALSE) && (status == TRUE))
							{
								status = FALSE;
							}
							else
							{
								status = FALSE;
							}
						}*/
			status &= st;
		}

		return status;
	}

	int 
	CActiveSection::ree_lower(TADVector& arg, TADVector& func, 
		TDoubleVector& err, TDoubleVector& out)
	{
		assert(arg.size() == err.size());

		int st = ree_lower(arg, err);

		out.clear();
		out.reserve(func.size());
		for (TADVector::iterator ifunc = func.begin(); ifunc != func.end(); ++ifunc)
		{
			assert(ifunc->pNode != NULL);
			out.push_back(ifunc->pNode->accum.d);
		}
		return st;
	}

} // namespace reed

#pragma warning(pop)	// restore state for each warning