void Assem_x86::assemLine( const string &line ) {

    int i=0;
    string name;
    vector<string> ops;

    //label?
    if( !isspace( line[i] ) ) {
        while( !isspace( line[i] ) ) ++i;
        string lab=line.substr( 0,i );
        if( !mod->addSymbol( lab.c_str(),mod->getPC() ) ) throw Ex( "duplicate label" );
    }

    //skip space
    while( isspace( line[i] ) && line[i]!='\n' ) ++i;
    if( line[i]=='\n' || line[i]==';' ) return;

    //fetch instruction name
    int from=i;
    for( ++i; !isspace( line[i] ); ++i ) {}
    name=line.substr( from,i-from );

    for(;;) {

        //skip space
        while( isspace( line[i] ) && line[i]!='\n' ) ++i;
        if( line[i]=='\n' || line[i]==';' ) break;

        int from=i;
        if( line[i]=='\"' ) {
            for( ++i; line[i]!='\"' && line[i]!='\n'; ++i ) {}
            if( line[i++]!='\"' ) throw Ex( "missing close quote" );
        } else {
            for( ++i; line[i]!=',' && line[i]!=';' && line[i]!='\n'; ++i ) {}
        }

        //back-up over space
        while( i && isspace( line[i-1] ) ) --i;
        ops.push_back( line.substr( from,i-from ) );

        //skip space
        while( isspace( line[i] ) && line[i]!='\n' ) ++i;
        if( line[i]=='\n' || line[i]==';' ) break;

        if( line[i++]!=',' ) throw Ex( "expecting ','" );
    }

    //pseudo op?
    if( name[0]=='.' ) {
        for( int k=0; k<ops.size(); ++k ) assemDir( name,ops[k] );
        return;
    }

    //normal instruction!
    if( ops.size()>2 ) throw Ex( "Too many operands" );
    ops.push_back( "" );
    ops.push_back( "" );
    assemInst( name,ops[0],ops[1] );
}

char* GFolderSerializer::nextPiece(size_t* pOutSize)
{
	switch(m_state)
	{
	case 0: // uncompressed header
		memcpy(m_pPos, "ugfs", 4);
		m_pPos += 4;
		m_size -= 4;
		m_state = 4;
		break;
	case 1: // figure out what to do next
		if(m_dirStack.size() > 0)
			continueDir();
		else
			return NULL;
		break;
	case 2: // continue reading file
		continueFile();
		break;
	case 3: // continue reading dir
		continueDir();
		break;
	case 4: // get started
		{
			// Change to the directory with the file or folder
			string sPath = m_szPath;
			if(sPath.length() > 0 && (sPath[sPath.length() - 1] == '/' || sPath[sPath.length() - 1] == '\\'))
				sPath.erase(sPath.length() - 1);
			PathData pd;
			GFile::parsePath(sPath.c_str(), &pd);
			if(pd.fileStart > 0)
			{
				string s;
				s.assign(m_szPath, pd.fileStart);
				if(chdir(s.c_str()) != 0)
					throw Ex("Failed to change dir to ", s.c_str());
			}

			// Add the file or folder
			if(access(m_szPath, 0) != 0)
				throw Ex("The file or folder ", m_szPath, " does not seem to exist");
			struct stat status;
			stat(m_szPath, &status);
			if(status.st_mode & S_IFDIR)
				startDir(m_szPath + pd.fileStart);
			else
				startFile(m_szPath + pd.fileStart);
		}
		break;
	default:
		throw Ex("Unexpected state");
	}
	*pOutSize = (m_pPos - m_pBuf);
	m_pPos = m_pBuf;
	m_size = BUF_SIZE;
	return m_pBuf;
}

	GPNGWriter()
	{
		m_pWriteStruct = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, error_handler, NULL);
		if(!m_pWriteStruct)
			throw Ex("Failed to create write struct. Out of mem?");
		m_pInfoStruct = png_create_info_struct(m_pWriteStruct);
		if(!m_pInfoStruct)
			throw Ex("Failed to create info struct. Out of mem?");
	}

// virtual
void GNaiveBayes::trainInner(const GMatrix& features, const GMatrix& labels)
{
	if(!features.relation().areNominal())
		throw Ex("GNaiveBayes only supports nominal features. Perhaps you should wrap it in a GAutoFilter.");
	if(!labels.relation().areNominal())
		throw Ex("GNaiveBayes only supports nominal labels. Perhaps you should wrap it in a GAutoFilter.");
	beginIncrementalLearningInner(features.relation(), labels.relation());
	for(size_t n = 0; n < features.rows(); n++)
		trainIncremental(features[n], labels[n]);
}

void FDTD2D::updateH(){
  // Bulk update, Hz
  for( int ii=0; ii<Nx-1; ii++){
    for( int jj=0; jj<Ny-1; jj++){
      Hz(ii,jj) += HzC(ii,jj)*(Ex(ii,jj+1) - Ex(ii,jj)
                             - Ey(ii+1,jj) + Ey(ii,jj));
    }
  }
  // No need to account for boundary conditions here; they are included in the E update equations.
  return;
}

// virtual
void GNaiveInstance::trainInner(const GMatrix& features, const GMatrix& labels)
{
	if(!features.relation().areContinuous())
		throw Ex("GNaiveInstance only supports continuous features. Perhaps you should wrap it in a GAutoFilter.");
	if(!labels.relation().areContinuous())
		throw Ex("GNaiveInstance only supports continuous labels. Perhaps you should wrap it in a GAutoFilter.");

	beginIncrementalLearningInner(features.relation(), labels.relation());
	for(size_t i = 0; i < features.rows(); i++)
		trainIncremental(features[i], labels[i]);
}

void GTokenizer::expect(const char* szString)
{
	while(*szString != '\0' && has_more())
	{
		char c = get();
		if(c != *szString)
			throw Ex("Expected \"", szString, "\" on line ", to_str(m_line), ", col ", to_str(col()));
		szString++;
	}
	if(*szString != '\0')
		throw Ex("Expected \", szString, \". Reached end-of-file instead.");
}

// virtual
void GLinearDistribution::trainInner(const GMatrix& features, const GMatrix& labels)
{
	if(!features.relation().areContinuous())
		throw Ex("GLinearDistribution only supports continuous features. Perhaps you should wrap it in a GAutoFilter.");
	if(!labels.relation().areContinuous())
		throw Ex("GLinearDistribution only supports continuous labels. Perhaps you should wrap it in a GAutoFilter.");

	// Init A with the inverse of the weights prior covariance matrix
	size_t dims = features.cols();
	GMatrix a(dims, dims);
	a.setAll(0.0);

	// Init XY
	size_t labelDims = labels.cols();
	GMatrix xy(dims, labelDims);
	xy.setAll(0.0);

	// Train on each instance
	double w = 1.0 / (m_noiseDev * m_noiseDev);
	for(size_t i = 0; i < features.rows(); i++)
	{
		// Update A
		const GVec& feat = features[i];
		for(size_t j = 0; j < dims; j++)
		{
			GVec& el = a[j];
			for(size_t k = 0; k < dims; k++)
				el[k] += feat[j] * feat[k];
		}

		// Update XY
		const GVec& lab = labels[i];
		for(size_t j = 0; j < dims; j++)
		{
			GVec& el = xy[j];
			for(size_t k = 0; k < labelDims; k++)
				el[k] += feat[j] * lab[k];
		}
	}
	a.multiply(w);
	xy.multiply(w);

	// Compute final matrices
	clear();
	m_pAInv = a.pseudoInverse();
	GAssert(m_pAInv->cols() == dims);
	GAssert(m_pAInv->rows() == dims);
	m_pWBar = GMatrix::multiply(xy, *m_pAInv, true, true);
	GAssert(m_pWBar->cols() == dims);
	GAssert(m_pWBar->rows() == labelDims);
	m_buf.resize(dims);
}

GDomNode* GKeyPair::serialize(GDom* pDoc, bool bIncludePrivateKey)
{
	GDomNode* pNode = pDoc->newObj();
	if(!n() || !publicKey())
		throw Ex("No key has been made yet");
	if(bIncludePrivateKey && !privateKey())
		throw Ex("This key-pair doesn't include the private key");
	pNode->addField(pDoc, "n", n()->serialize(pDoc));
	pNode->addField(pDoc, "public", publicKey()->serialize(pDoc));
	if(bIncludePrivateKey)
		pNode->addField(pDoc, "private", privateKey()->serialize(pDoc));
	return pNode;
}

void GNaiveBayes_CheckResults(double yprior, double ycond, double nprior, double ncond, GPrediction* out)
{
	double py = yprior * ycond;
	double pn = nprior * ncond;
	double sum = py + pn;
	py /= sum;
	pn /= sum;
	GCategoricalDistribution* pCat = out->asCategorical();
	GVec& vals = pCat->values(2);
	if(std::abs(vals[0] - py) > 1e-8)
		throw Ex("wrong");
	if(std::abs(vals[1] - pn) > 1e-8)
		throw Ex("wrong");
}

bool SearchCProcs(var &result, var head, var body)
{
	stdext::hash_map<Var,CProc>::const_iterator
		iter = CProcs.find(head);
	if(iter != CProcs.end())
	{
		PushCProcSymbol pcs(head);
		if (Verbose >= 3) {
			wcerr << _W("Enter cproc for ");
			Println(Ex(head, body), wcerr);
		}
		try {
			result = iter->second(body);
		}
		catch (UserException&) {
			throw;
		}
// 		catch (std::exception& e) {
// 			wcerr << L"std::exception thrown: " << mbs2wcs(e.what()) << std::endl;
// 			wcerr << "CProc error while eval ";
// 			Println(Ex(head, body), wcerr);
// 			return false;
// 		}
		catch (...) {
			wcerr << "CProc error while eval ";
			Println(Ex(head, body), wcerr);
			throw;
		}
		if (result)
		{
			if (TraceRuleSymbols.count(head) != 0)
			{
				// FIXME: rewrite this after we implement the kernel message infrastructure
				Print(head);
				wcout << _W("::tracer : ");
				Print(Ex(head, body));
				wcout << _W(" :> ");
				Println(result);
			}
			if (Verbose >= 3) {
				wcerr << _W("Result is ");
				Println(result, wcerr);
			}
			return true;
		}
	}
	return false;
}

void GImage::loadPng(const char* szFilename)
{
	size_t nSize;
	std::ifstream s;
	char* pBuf;
	s.exceptions(std::ios::badbit);
	try
	{
		// NB: we double copy the content of the file!
		// Once in a std::string and then to final destination
		// We do this so it will work with named pipes, which
		// do not permit seeking to the end of the file to determine
		// the file length.
		s.open(szFilename, std::ios::binary);
		std::stringstream buffer;
		buffer << s.rdbuf();
		std::string content = buffer.str();
		nSize = content.size();
		pBuf = new char[nSize];
		memcpy(pBuf, content.data(), nSize);
		s.close();
	}
	catch(const std::exception&)
	{
		throw Ex("Error while trying to open the file, ", szFilename, ". ", strerror(errno));
	}
	readPng(this, (const unsigned char*)pBuf, nSize);
	delete[] pBuf;
}

void GFolderSerializer::continueDir()
{
	GDirList* pDL = m_dirStack.top();
	if(pDL->m_folders.size() > 0)
	{
		startDir(pDL->m_folders.back().c_str());
		pDL->m_folders.pop_back();
	}
	else if(pDL->m_files.size() > 0)
	{
		startFile(pDL->m_files.back().c_str());
		pDL->m_files.pop_back();
	}
	else
	{
		// End of dir indicator
		*m_pPos = 'e';
		m_pPos++;
		m_size--;

		// Move out of the dir
		delete(pDL);
		m_dirStack.pop();
		if(chdir("..") != 0)
			throw Ex("Failed to chdir to ..");
		m_state = 1;
	}
}

void GFolderSerializer::startFile(const char* szFilename)
{
	// File indicator
	*m_pPos = 'f';
	m_pPos++;
	m_size--;
	addName(szFilename);

	// The file size
	m_pInStream = new std::ifstream();
	unsigned long long size = 0;
	try
	{
		m_pInStream->exceptions(std::ios::badbit | std::ios::failbit);
		m_pInStream->open(szFilename, std::ios::binary);
		m_pInStream->seekg(0, std::ios::end);
		size = m_pInStream->tellg();
		m_pInStream->seekg(0, std::ios::beg);
	}
	catch(const std::exception e)
	{
		throw Ex("Error opening file: ", szFilename);
	}
	memcpy(m_pPos, &size, sizeof(unsigned long long));
	m_pPos += sizeof(unsigned long long);
	m_size -= sizeof(unsigned long long);
	m_remaining = (size_t)size;

	// The file
	continueFile();
}

GPlotLabelSpacer::GPlotLabelSpacer(double min, double max, int maxLabels)
{
	if(max <= min)
		throw Ex("invalid range");
	int p = (int)ceil(log((max - min) / maxLabels) * M_LOG10E);

	// Every 10
	m_spacing = pow(10.0, p);
	m_start = (int)ceil(min / m_spacing);
	m_count = (int)floor(max / m_spacing) - m_start + 1;

	if(m_count * 5 + 4 < maxLabels)
	{
		// Every 2
		m_spacing *= 0.2;
		m_start = (int)ceil(min / m_spacing);
		m_count = (int)floor(max / m_spacing) - m_start + 1;
	}
	else if(m_count * 2 + 1 < maxLabels)
	{
		// Every 5
		m_spacing *= 0.5;
		m_start = (int)ceil(min / m_spacing);
		m_count = (int)floor(max / m_spacing) - m_start + 1;
	}
}

Type objective_function<Type>::operator() ()
{
	DATA_VECTOR(times);
	DATA_VECTOR(obs);
	
	PARAMETER(log_R0);
	PARAMETER(m);
	PARAMETER(log_theta);
	PARAMETER(log_sigma);
	Type theta=exp(log_theta);
	Type sigma=exp(log_sigma);
	Type R0=exp(log_R0);
	
	int n1=times.size();
	int n2=2;//mean and variance
	vector<Type> Dt(n1-1);
	vector<Type> Ex(n1-1);
	vector<Type> Vx(n1-1);
	Type nll=0;
	m=0;
	Dt=diff(times);
	Ex=theta*(Type(1)-exp(-R0*Dt)) + obs.segment(0, n1-1)*exp(-R0*Dt);
	Vx=Type(0.5)*sigma*sigma*(Type(1)-exp(Type(-2)*R0*Dt))/R0;
	
	for(int i=0; i<n1-1; i++)
	{
		nll-= dnorm(obs[i+1], Ex[i], sqrt(Vx[i]), true);
	}
	return nll;
}

void Assem_x86::emitImm( const string &s,int size ) {

    Operand op(s);
    op.parse();
    if( !(op.mode&IMM) ) throw Ex( "operand must be immediate" );
    emitImm( op,size );
}

void CommReg::SendFields(UInt nfields, MEField<> *const *sfields, MEField<> *const *rfields) {
  // Get all the subfields and send out over the specs.
  std::vector<_field*> sf;
  std::vector<_field*> rf;
  UInt obj_type = 0;
  for (UInt i = 0; i < nfields; i++) {
    if (&sfields[i]->GetMEFamily() != &rfields[i]->GetMEFamily())
      throw Ex() << "Send fields, me for " << sfields[i]->name() << " does not match rfield:"
                 << rfields[i]->name();
    sfields[i]->Getfields(sf);
    rfields[i]->Getfields(rf);
  }

  for (UInt j = 0; j < sf.size(); j++)
    obj_type |= sf[j]->GetAttr().GetType();

/*
std::cout << "sf size=" << sf.size() << ". _fields are:";
for (UInt i = 0; i < sf.size(); i++) {
  std::cout << sf[i]->name() << ", dim=" << sf[i]->dim() << std::endl;
}
std::cout << "rf size=" << rf.size() << ". _fields are:";
for (UInt i = 0; i < rf.size(); i++) {
  std::cout << rf[i]->name() << ", dim=" << rf[i]->dim() << std::endl;
}
*/
  ThrowRequire(sf.size() == rf.size());
  // Now send via the spec(s)
  // TODO: be smarter: select only the relevant spec to send each field.
  if ((obj_type & MeshObj::NODE)) node_rel.send_fields(sf.size(), &sf[0], &rf[0]);
  if ((obj_type & MeshObj::EDGE)) edge_rel.send_fields(sf.size(), &sf[0], &rf[0]);
  if ((obj_type & MeshObj::FACE)) face_rel.send_fields(sf.size(), &sf[0], &rf[0]);
  if ((obj_type & MeshObj::ELEMENT)) elem_rel.send_fields(sf.size(), &sf[0], &rf[0]);
}

// static
GActivationFunction* GActivationFunction::deserialize(GDomNode* pNode)
{
	const char* szName = pNode->asString();
	if(strcmp(szName, "logistic") == 0)
		return new GActivationLogistic();
	else if(strcmp(szName, "arctan") == 0)
		return new GActivationArcTan();
	else if(strcmp(szName, "tanh") == 0)
		return new GActivationTanH();
	else if(strcmp(szName, "algebraic") == 0)
		return new GActivationAlgebraic();
	else if(strcmp(szName, "identity") == 0)
		return new GActivationIdentity();
	else if(strcmp(szName, "gaussian") == 0)
		return new GActivationGaussian();
	else if(strcmp(szName, "bidir") == 0)
		return new GActivationBiDir();
	else if(strcmp(szName, "bend") == 0)
		return new GActivationBend();
	else if(strcmp(szName, "sinc") == 0)
		return new GActivationSinc();
	else if(strcmp(szName, "piecewise") == 0)
		return new GActivationPiecewise();
	else
		throw Ex("Unrecognized activation function: ", szName);
	return NULL;
}

/// Refines this model based on a recently performed action and change in beliefs
void TransitionModel::trainIncremental(const GVec& beliefs, const GVec& actions, const GVec& nextBeliefs)
{
	// Buffer the pattern
	GVec& destIn = trainInput.row(trainPos);
	GVec& destOut = trainOutput.row(trainPos);
	trainPos++;
	trainSize = std::max(trainSize, trainPos);
	if(trainPos >= trainInput.rows())
		trainPos = 0;
	if(beliefs.size() + actions.size() != destIn.size() || beliefs.size() != destOut.size())
		throw Ex("size mismatch");
	destIn.put(0, beliefs);
	destIn.put(beliefs.size(), actions);
	for(size_t i = 0; i < destOut.size(); i++)
		destOut[i] = 0.5 * (nextBeliefs[i] - beliefs[i]);
/*
destIn.print();
std::cout << "->";
destOut.print();
std::cout << "\n";
std::cout << to_str(0.5 * cos(destIn[2])) << ", " << to_str(0.5 * sin(destIn[2])) << "\n";
*/
	// Refine the model
	size_t iters = std::min(trainIters, 1000 * trainSize);
	for(size_t i = 0; i < iters; i++)
		doSomeTraining();
}