/*
 * \file DataSet.c
 * \brief: implementation of the DataSet class
 */


#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <vector>
#include <functional>
#include <algorithm>
#include <iostream>

#include "./DataSet.h"
#include "../shared/shared.h"
#include "../include/macros.h"
#include "../EnumDefinitions/EnumDefinitions.h"

using namespace std;

DataSet::DataSet(){
	
	sorted=0;
	sorted_ids=NULL;
	id_offsets=NULL;

}

DataSet::DataSet(int dataset_enum){
	enum_type=dataset_enum;
	
	sorted=0;
	sorted_ids=NULL;
	id_offsets=NULL;

}

DataSet::~DataSet(){
	clear();
	xfree((void**)&sorted_ids);
	xfree((void**)&id_offsets);
}

int  DataSet::GetEnum(){
	return enum_type;
}
	
#undef __FUNCT__
#define  __FUNCT__ "DataSet::Echo"

void DataSet::Echo(){

	
	vector<Object*>::iterator object;

	if(this==NULL)throw ErrorException(__FUNCT__," trying to echo a NULL dataset");

	_printf_("DataSet echo: %i objects\n",objects.size());

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Call echo on object: */
		(*object)->Echo();

	}
	return;
}

int DataSet::MarshallSize(){

	vector<Object*>::iterator object;
	int                      marshalled_dataset_size=0;


	for ( object=objects.begin() ; object < objects.end(); object++ ){
		marshalled_dataset_size+= (*object)->MarshallSize();
	}

	marshalled_dataset_size+=sizeof(int); //objects size
	marshalled_dataset_size+=sizeof(int); //sorted size
	if(sorted){
		marshalled_dataset_size+=(int)objects.size()*sizeof(int); //sorted ids
		marshalled_dataset_size+=(int)objects.size()*sizeof(int); //id offsets
	}

	return marshalled_dataset_size;
}
#undef __FUNCT__
#define __FUNCT__ "DataSet::Marshall"

char* DataSet::Marshall(){

	vector<Object*>::iterator object;
	int                       object_size;
	int                       marshalled_dataset_size=0;
	char*                     marshalled_dataset=NULL;
	char*                     old_marshalled_dataset=NULL;

	/*First get size of marshalled dataset: */
	object_size=(int)objects.size();

	marshalled_dataset_size=MarshallSize();
	
	/*Allocate marshalled dataset: */
	marshalled_dataset=(char*)xmalloc(marshalled_dataset_size*sizeof(char)); 

	/*Keep track of old_marshalled_dataset: */
	old_marshalled_dataset=marshalled_dataset;

	/*Store internals of dataset first: */
	memcpy(marshalled_dataset,&object_size,sizeof(int)); marshalled_dataset+=sizeof(int);
	memcpy(marshalled_dataset,&sorted,sizeof(int)); marshalled_dataset+=sizeof(int);
	if(sorted){
		if(object_size)memcpy(marshalled_dataset,sorted_ids,object_size*sizeof(int)); marshalled_dataset+=object_size*sizeof(int);
		if(object_size)memcpy(marshalled_dataset,id_offsets,object_size*sizeof(int)); marshalled_dataset+=object_size*sizeof(int);
	}

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		(*object)->Marshall(&marshalled_dataset);
	}

	/* Ok, marshalled_dataset now points to the end of the original marshalled_dataset pointer 
	 * before  we started the loop on objects. Get object to point right again: */
	marshalled_dataset-=marshalled_dataset_size;

	/*We should be back to old_marshalled_dataset: check and abort if that's not the case, 
	 * because this is a nasty error: */
	if (marshalled_dataset!=old_marshalled_dataset){
		_printf_("%s%s\n",__FUNCT__," error message: final marshalled dataset is different from initial one!");
		abort();
	}

	/*Return: */
	return marshalled_dataset;
}

int  DataSet::AddObject(Object* object){

	objects.push_back(object);
	
	return 1;
}

int  DataSet::DeleteObject(int id){

	return 0;
}

int  DataSet::Size(void){

	return objects.size();
}

#undef __FUNCT__
#define __FUNCT__ "DataSetDemarshall"
DataSet* DataSetDemarshall(char* marshalled_dataset){

	int i;
	
	DataSet* dataset=NULL;
	int      numobjects=0;
	int      enum_type;
	Object*  object=NULL;
	int      sorted;
	int*     sorted_ids=NULL;
	int*     id_offsets=NULL;

	/*initialize dataset: */
	dataset=new DataSet();

	/*Get internals first: */
	memcpy(&numobjects,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
	memcpy(&sorted,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
	if(sorted){
		if(numobjects){
			sorted_ids=(int*)xmalloc(numobjects*sizeof(int));
			id_offsets=(int*)xmalloc(numobjects*sizeof(int));
			memcpy(sorted_ids,marshalled_dataset,numobjects*sizeof(int)); marshalled_dataset+=numobjects*sizeof(int);
			memcpy(id_offsets,marshalled_dataset,numobjects*sizeof(int)); marshalled_dataset+=numobjects*sizeof(int);
		}
		dataset->SetSorting(sorted_ids,id_offsets);
	}

	#ifdef _DEBUG_
	_printf_("Number of objects in dataset being demarshalled: %i\n",numobjects);
	#endif

	for(i=0;i<numobjects;i++){
	
		/*get enum type of object: */
		memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);

		if(enum_type==NodeEnum()){
			Node* node=NULL;
			node=new Node();
			node->Demarshall(&marshalled_dataset);
			dataset->AddObject(node);
		}
		else if(enum_type==ParamEnum()){
			Param* param=NULL;
			param=new Param();
			param->Demarshall(&marshalled_dataset);
			dataset->AddObject(param);
		}
		else if(enum_type==TriaEnum()){
			Tria* tria=NULL;
			tria=new Tria();
			tria->Demarshall(&marshalled_dataset);
			dataset->AddObject(tria);
		}
		else if(enum_type==PentaEnum()){
			Penta* penta=NULL;
			penta=new Penta();
			penta->Demarshall(&marshalled_dataset);
			dataset->AddObject(penta);
		}
		else if(enum_type==MaticeEnum()){
			Matice* matice=NULL;
			matice=new Matice();
			matice->Demarshall(&marshalled_dataset);
			dataset->AddObject(matice);
		}
		else if(enum_type==MatparEnum()){
			Matpar* matpar=NULL;
			matpar=new Matpar();
			matpar->Demarshall(&marshalled_dataset);
			dataset->AddObject(matpar);
		}
		else if(enum_type==SpcEnum()){
			Spc* spc=NULL;
			spc=new Spc();
			spc->Demarshall(&marshalled_dataset);
			dataset->AddObject(spc);
		}
		else if(enum_type==PenpairEnum()){
			Penpair* penpair=NULL;
			penpair=new Penpair();
			penpair->Demarshall(&marshalled_dataset);
			dataset->AddObject(penpair);
		}
		else if(enum_type==IcefrontEnum()){
			Icefront* icefront=NULL;
			icefront=new Icefront();
			icefront->Demarshall(&marshalled_dataset);
			dataset->AddObject(icefront);
		}
		else if(enum_type==RgbEnum()){
			Rgb* rgb=NULL;
			rgb=new Rgb();
			rgb->Demarshall(&marshalled_dataset);
			dataset->AddObject(rgb);
		}
		else{
			throw ErrorException(__FUNCT__,exprintf("%s%i"," could not recognize enum type: ",enum_type));
		}

	}
	return dataset;

}
		
int   DataSet::FindParam(void* pvalue, char* name){

	/*Go through a dataset, and find a Param* object 
	 *which parameter name is "name" : */
	
	vector<Object*>::iterator object;
	Param* param=NULL;

	int found=0;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Find param type objects: */
		if((*object)->Enum()==ParamEnum()){

			/*Ok, this object is a parameter, recover it and ask which name it has: */
			param=(Param*)(*object);

			if (strcmp(param->GetParameterName(),name)==0){
				/*Ok, this is the one! Recover the value of this parameter: */
				param->GetParameterValue(pvalue);
				found=1;
				break;
			}
		}
	}
	return found;
}

void   DataSet::NodeRank(int* ranks){

	/*Go through a dataset, and find a Node* object. For this object, 
	 * ask it to report its cpu rank: */
	
	int node_rank;
	int id;
	vector<Object*>::iterator object;
	
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Call echo on object: */
		if((*object)->Enum()==NodeEnum()){
			
			/*Ok, this object is a node, ask which rank it has: */
			node_rank=(*object)->MyRank();
			
			/*Which id does it have: */
			id=(*object)->GetId();

			/*Plug rank into ranks, according to id: */
			ranks[id-1]=node_rank; //Matlab ids start at 1.
		}
	}
	return;
}


void  DataSet::DistributeDofs(int numberofnodes,int numdofspernode){

	int  dofcount=0;
	int* alldofcount=NULL;
	int* borderdofs=NULL;
	int* allborderdofs=NULL;
	int  i;
	vector<Object*>::iterator object;
	Node* node=NULL;

	extern int num_procs;
	extern int my_rank;
	
	dofcount=0;
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){

			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to distribute dofs, and update the dofcount: */
			node->DistributeDofs(&dofcount);
			
		}
	}

	/*Ok, now every node has distributed dofs, but locally, and with a dof count starting from 
	 *0. This means the dofs between all the cpus are not synchronized! We need to synchronize all 
	 *dof on all cpus, and use those to update the dofs of every node: */
	
	alldofcount=(int*)xmalloc(num_procs*sizeof(int));
	MPI_Gather(&dofcount,1,MPI_INT,alldofcount,1,MPI_INT,0,MPI_COMM_WORLD);
	MPI_Bcast(alldofcount,num_procs,MPI_INT,0,MPI_COMM_WORLD);

	/*Ok, now every cpu should start its own dof count at the end of the dofcount 
	 * from cpu-1. : */
	dofcount=0;
	if(my_rank==0){
		dofcount=0;
	}
	else{
		for(i=0;i<my_rank;i++){
			dofcount+=alldofcount[i];
		}
	}


	/*Ok, now every cpu knows where his dofs should start. Update the dof count: */
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){

			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to update his dofs: */
			node->UpdateDofs(dofcount);
			
		}
	}

	/*Finally, remember that cpus may have skipped some nodes, when these nodes were 
	 * clones: */
	borderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
	allborderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){

			node=(Node*)(*object);
			
			/*Ok, let this object show its border dofs, if is is a border dof: */
			node->ShowBorderDofs(borderdofs);
			
		}
	}
	MPI_Allreduce ( (void*)borderdofs,(void*)allborderdofs,numberofnodes*numdofspernode,MPI_INT,MPI_MAX,MPI_COMM_WORLD);

	/*Ok, now every cpu knows everyone else's border node dofs, update the border dofs accordingly: */
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){

			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to update his dofs: */
			node->UpdateBorderDofs(allborderdofs);
			
		}
	}


	/* Free ressources: */
	xfree((void**)&alldofcount);
	xfree((void**)&borderdofs);
	xfree((void**)&allborderdofs);
	return;

}

		
void  DataSet::CreatePartitioningVector(Vec* ppartition,int numberofnodes,int numdofspernode){

	/*output: */
	Vec partition=NULL;
	vector<Object*>::iterator object;
	Node* node=NULL;

	/*Create partition vector: */
	partition=NewVec(numberofnodes*numdofspernode);

	/*Go through all nodes, and ask each node to plug its doflist in 
	 * partition. The location where each node plugs its doflist into 
	 * partition is determined by its (id-1)*numdofspernode (ie, serial 
	 * organisation of the dofs).
	 */
	
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){
			
			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to plug values into partition: */
			node->CreatePartition(partition);
			
		}
	}

	/*Assemble the petsc vector: */
	VecAssemblyBegin(partition);
	VecAssemblyEnd(partition);

	/*Assign output pointers: */
	*ppartition=partition;

	return;
}

void  DataSet::FlagClones(int numberofnodes){

	int i;
	extern int num_procs;
	
	int* ranks=NULL;
	int* minranks=NULL;

	vector<Object*>::iterator object;
	Node* node=NULL;
	
	/*Allocate ranks: */
	ranks=(int*)xmalloc(numberofnodes*sizeof(int));
	minranks=(int*)xmalloc(numberofnodes*sizeof(int));
	
	for(i=0;i<numberofnodes;i++)ranks[i]=num_procs; //no cpu can have rank num_procs. This is the maximum limit.

	/*Now go through all our nodes and ask them to report to who they belong (which rank): */
	NodeRank(ranks);

	#ifdef _DEBUG_
	for(i=0;i<numberofnodes;i++){
		_printf_("%i\n",ranks[i]);
	}
	#endif

	/*We need to take the minimum rank for each node, and every cpu needs to get that result. That way, 
	 * when we start building the dof list for all nodes, a cpu can check whether its node already has been 
	 * dealt with by another cpu. We take the minimum because we are going to manage dof assignment in increasing 
	 * order of cpu rank. This is also why we initialized this array to num_procs.*/
	MPI_Allreduce ( (void*)ranks,(void*)minranks,numberofnodes,MPI_INT,MPI_MIN,MPI_COMM_WORLD);

	#ifdef _DEBUG_
	for(i=0;i<numberofnodes;i++){
		_printf_("%i\n",minranks[i]);
	}
	#endif

	/*Now go through all nodes, and use minranks to flag which nodes are cloned: */
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){
			
			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to plug values into partition: */
			node->SetClone(minranks);
			
		}
	}

	/*Free ressources: */
	xfree((void**)&ranks); 
	xfree((void**)&minranks);

}
		

int   DataSet::NumberOfDofs(){


	vector<Object*>::iterator object;
	Node* node=NULL;
	int   numdofs=0;
	int   allnumdofs;

	/*Now go through all nodes, and get how many dofs they own, unless they are clone nodes: */
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a node: */
		if((*object)->Enum()==NodeEnum()){
			
			node=(Node*)(*object);
			
			/*Ok, this object is a node, ask it to plug values into partition: */
			if (!node->IsClone()){

				numdofs+=node->GetNumberOfDofs();
				
			}
		}
	}

	/*Gather from all cpus: */
	MPI_Allreduce ( (void*)&numdofs,(void*)&allnumdofs,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

	return allnumdofs;
}

#undef __FUNCT__ 
#define __FUNCT__ "DataSet::SetupSpcs"
		
void   DataSet::SetupSpcs(DataSet* nodes,Vec yg){

	vector<Object*>::iterator object;
	Spc* spc=NULL;
	Node* node=NULL;

	int nodeid;
	int dof;
	double value;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a single point constraint (spc): */
		if((*object)->Enum()==SpcEnum()){
			
			spc=(Spc*)(*object);
			
			/*Ok, this object is a constraint. Get the nodeid from the node it applies to: */
			nodeid=spc->GetNodeId();
			dof=spc->GetDof();
			value=spc->GetValue();

			/*Now, chase through nodes and find the corect node: */
			node=(Node*)nodes->GetObjectById(NULL,nodeid);

			/*Apply constraint: */
			if(node){ //in case the spc is dealing with a node on another cpu
				node->ApplyConstraint(yg,dof,value);
			}

		}
	}

	/*Assemble yg: */
	VecAssemblyBegin(yg);
	VecAssemblyEnd(yg);
}
		

		
int   DataSet::NumberOfRgbs(){

	vector<Object*>::iterator object;
	Rgb* rgb=NULL;
	int  count=0;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a single point constraint (spc): */
		if((*object)->Enum()==RgbEnum()){ //we assume uniqueness of all Rgbs, no error checking here.

			count++;

		}
	}

	return count;
}

#undef __FUNCT__ 
#define __FUNCT__  "DataSet::SetupMpcs"
void DataSet::SetupMpcs(Mat Rmg,DataSet* nodes){

	vector<Object*>::iterator object;
	Rgb* rgb=NULL;
	int  nodeid1;
	int  nodeid2;
	int  dof;

	int  dof1;
	int  dof2;


	Node* node1=NULL;
	Node* node2=NULL;

	int count=-1;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a single point constraint (spc): */
		if((*object)->Enum()==RgbEnum()){ //we assume uniqueness of all Rgbs, no error checking here.

			/*we found an rgb, increment counter, so that row index for Rmg is up to date: */
			count++;
			rgb=(Rgb*)(*object);


			nodeid1=rgb->GetNodeId1();
			nodeid2=rgb->GetNodeId2();
			dof=rgb->GetDof();

			/*For this rgb, find the nodes that go with it: */
			node1=(Node*)nodes->GetObjectById(NULL,nodeid1);
			node2=(Node*)nodes->GetObjectById(NULL,nodeid2);

			if ((node1 && !node2) || (!node1 && node2)){
				/*we are missing one node, not good!*/
				throw ErrorException(__FUNCT__,exprintf("%s%p%s%p"," in Rgb, missing one node. node1: ",node1," node2: ",node2));
			}

			if(!node1 && !node2){
				/*That's ok, this Rgb can't find those nodes, so leave them alone. They are probably not on this 
				 * cpu!*/
			}
			else{
				/*Ok, this cpu owns both nodes. Put dof for node1 into m set: */
				node1->DofInMSet(dof-1);
				

				/*Plug values into Rmg. We essentially want dofs from node1 and node2 to be the 
				 *same: */
				dof1=node1->GetDof(dof-1); //matlab indexing
				dof2=node2->GetDof(dof-1); //matlab indexing

				MatSetValue(Rmg,count,dof1,1.0,INSERT_VALUES);
				MatSetValue(Rmg,count,dof2,-1.0,INSERT_VALUES);

			}
		}
	}
}
	
void DataSet::FlagNodeSets(Vec pv_g, Vec pv_m, Vec pv_n, Vec pv_f, Vec pv_s){

	vector<Object*>::iterator object;
	Node* node=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Check this is a single point constraint (spc): */
		if((*object)->Enum()==NodeEnum()){
			
			node=(Node*)(*object);

			/*Plug set values intp our 4 set vectors: */
			node->CreateVecSets(pv_g,pv_m,pv_n,pv_f,pv_s);

		}
	}

	/*Assemble: */
	VecAssemblyBegin(pv_g);
	VecAssemblyEnd(pv_g);

	VecAssemblyBegin(pv_m);
	VecAssemblyEnd(pv_m);

	VecAssemblyBegin(pv_n);
	VecAssemblyEnd(pv_n);

	VecAssemblyBegin(pv_f);
	VecAssemblyEnd(pv_f);

	VecAssemblyBegin(pv_s);
	VecAssemblyEnd(pv_s);

}
void  DataSet::clear(){
	
	vector<Object*>::iterator object;
	
	for ( object=objects.begin() ; object < objects.end(); object++ ){
		delete (*object);
		object=objects.erase(object);
	}
}


void DataSet::Configure(DataSet* elements,DataSet* loads, DataSet* nodes, DataSet* materials){

	vector<Object*>::iterator object;
	Element* element=NULL;
	Load* load=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsElement((*object)->Enum())){


			element=(Element*)(*object);

			element->Configure(loads,nodes,materials);
		}
		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);

			load->Configure(elements,nodes,materials);
		}

	}

}

void DataSet::Presort(){
	
	/*vector of objects is already sorted, just allocate the sorted ids and their
	 * offsets:*/
	int i;

	if(objects.size()){
		sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
		id_offsets=(int*)xmalloc(objects.size()*sizeof(int));
		for(i=0;i<objects.size();i++){
			id_offsets[i]=i;
			sorted_ids[i]=objects[i]->GetId();
		}
	}

	/*set sorted flag: */
	sorted=1;
}

#undef __FUNCT__
#define __FUNCT__ "DataSet::Sort"
void DataSet::Sort(){

	/*Only sort if we are not already sorted: */
	if(!sorted){
		throw ErrorException(__FUNCT__," not implemented yet!");
	}
}


Object* DataSet::GetObjectByOffset(int offset){

	return objects[offset];

}



#undef __FUNCT__ 
#define __FUNCT__ "GetObjectById"
Object* DataSet::GetObjectById(int* poffset,int eid){

	int id_offset;
	int offset;
	int i;

	if(!sorted)throw ErrorException(__FUNCT__," trying to binary search on a non-sorted dataset!");

	/*Carry out a binary search on the sorted_ids: */
	if(!binary_search(&id_offset,eid, sorted_ids,objects.size())){
		throw ErrorException(__FUNCT__,exprintf("%s%i"," could not find object with id ",eid));
	}

	/*Convert  the id offset into sorted offset: */
	offset=id_offsets[id_offset];

	/*Assign output pointers if requested:*/
	if (poffset)*poffset=offset;

	/*Return object at offset position in objects :*/
	return objects[offset];
}


#undef __FUNCT__
#define __FUNCT__ "DataSet::SetSorting"

void DataSet::SetSorting(int* in_sorted_ids,int* in_id_offsets){

	sorted=1;
	sorted_ids=in_sorted_ids;
	id_offsets=in_id_offsets;
}

#undef __FUNCT__
#define __FUNCT__ "DataSet::CreateKMatrix"
void  DataSet::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){

	vector<Object*>::iterator object;
	Element* element=NULL;
	Load* load=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			element->CreateKMatrix(Kgg,inputs,analysis_type);
		}
		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);
			load->CreateKMatrix(Kgg,inputs,analysis_type);
		}
	}

}

#undef __FUNCT__
#define __FUNCT__ "DataSet::CreatePVector"
void  DataSet::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){

	vector<Object*>::iterator object;
	Element* element=NULL;
	Load* load=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			element->CreatePVector(pg,inputs,analysis_type);
		}
		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);
			load->CreatePVector(pg,inputs,analysis_type);
		}		
	}

}

#undef __FUNCT__		
#define __FUNCT__ "UpdateFromInputs"
void  DataSet::UpdateFromInputs(ParameterInputs* inputs){

	vector<Object*>::iterator object;
	Element* element=NULL;
	Load* load=NULL;
	Material* material=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			element->UpdateFromInputs(inputs);
		}
		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);
			load->UpdateFromInputs(inputs);
		}
		if(EnumIsMaterial((*object)->Enum())){

			material=(Material*)(*object);
			material->UpdateFromInputs(inputs);
		}		
	}

}


#undef __FUNCT__
#define __FUNCT__ "DataSet::PenaltyCreateKMatrix"
void  DataSet::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){

	vector<Object*>::iterator object;
	Load* load=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);
			load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type);
		}
	}

}

#undef __FUNCT__
#define __FUNCT__ "DataSet::PenaltyCreatePVector"
void  DataSet::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){

	vector<Object*>::iterator object;
	Load* load=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){

		if(EnumIsLoad((*object)->Enum())){

			load=(Load*)(*object);
			load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type);
		}		
	}

}

#undef __FUNCT__
#define __FUNCT__ "RiftIsPresent"
int   DataSet::RiftIsPresent(){

	return 0; //for now

}

#undef __FUNCT__
#define __FUNCT__ "RiftConstraints"
void  DataSet::RiftConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){
	
	throw ErrorException(__FUNCT__," not implemented yet!");

}

#undef __FUNCT__
#define __FUNCT__ "MeltingIsPresent"
int   DataSet::MeltingIsPresent(){


	return 0; //for now
	
}

#undef __FUNCT__
#define __FUNCT__ "MeltingConstraints"
void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){

	throw ErrorException(__FUNCT__," not implemented yet!");

}

#undef __FUNCT__
#define __FUNCT__ "DataSet::Copy"
DataSet*   DataSet::Copy(void){


	DataSet* copy=NULL;
	vector<Object*>::iterator object;
	Object* object_copy=NULL;

	copy=new DataSet(enum_type);

	copy->sorted=sorted;
	copy->presorted=presorted;
	if(sorted_ids){
		copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
		memcpy(copy->sorted_ids,sorted_ids,objects.size()*sizeof(int));
	}
	if(id_offsets){
		copy->sorted_ids=(int*)xmalloc(objects.size()*sizeof(int));
		memcpy(copy->sorted_ids,id_offsets,objects.size()*sizeof(int));
	}

	/*Now we need to deep copy the objects: */
	for ( object=objects.begin() ; object < objects.end(); object++ ){

		/*Call echo on object: */
		object_copy = (*object)->copy();
		copy->AddObject(object_copy);
	}
	return copy;
}


void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){


	vector<Object*>::iterator object;
	Element* element=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			element->Du(du_g,u_g,u_g_obs,inputs,analysis_type);
		}
	}


}

void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type){


	vector<Object*>::iterator object;
	Element* element=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,control_type);
		}
	}


}		
		
void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){

	double J=0;;
	
	vector<Object*>::iterator object;
	Element* element=NULL;

	for ( object=objects.begin() ; object < objects.end(); object++ ){
		
		if(EnumIsElement((*object)->Enum())){

			element=(Element*)(*object);
			J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type);

		}
	}

	/*Assign output pointers:*/
	*pJ=J;

}
