diff --git a/BlockTriangular.cc b/BlockTriangular.cc index b9348f6d669cedb6f2b2d051d59444b3f31e3550..5f9b02c3f76367032304fad532f82f87a69bbf0e 100644 --- a/BlockTriangular.cc +++ b/BlockTriangular.cc @@ -108,19 +108,18 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block bool *tmp_variable_evaluated; bool *Cur_IM; bool *IM, OK; - ModelBlock->Periods = periods; int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo; + ModelBlock->Periods = periods; ModelBlock->Block_List[count_Block].is_linear=true; ModelBlock->Block_List[count_Block].Size = size; ModelBlock->Block_List[count_Block].Type = type; - ModelBlock->Block_List[count_Block].Temporary_terms=new temporary_terms_type (); - ModelBlock->Block_List[count_Block].Temporary_terms->clear(); ModelBlock->Block_List[count_Block].Temporary_InUse=new temporary_terms_inuse_type (); ModelBlock->Block_List[count_Block].Temporary_InUse->clear(); ModelBlock->Block_List[count_Block].Simulation_Type = SimType; ModelBlock->Block_List[count_Block].Equation = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); ModelBlock->Block_List[count_Block].Variable = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); + ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation = (temporary_terms_type**)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(temporary_terms_type)); ModelBlock->Block_List[count_Block].Own_Derivative = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int)); Lead = Lag = 0; first_count_equ = *count_Equ; @@ -144,6 +143,8 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool)); for (i = 0;i < size;i++) { + ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]=new temporary_terms_type (); + ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]->clear(); ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[count_Block].Variable[i] = Index_Var_IM[*count_Equ].index; i_1 = Index_Var_IM[*count_Equ].index; @@ -453,7 +454,9 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const } } free(ModelBlock->Block_List[blk].IM_lead_lag); - delete(ModelBlock->Block_List[blk].Temporary_terms); + for(i=0; i<ModelBlock->Block_List[blk].Size; i++) + delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i]; + free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation); delete(ModelBlock->Block_List[blk].Temporary_InUse); } free(ModelBlock->Block_List); diff --git a/ExprNode.cc b/ExprNode.cc index e9433b37a21c7a36f670fd5e57d186d4bdb8bc98..79b1b9ac570983ad46d8959fe10828affa759016 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -84,9 +84,10 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count, void ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const { // Nothing to do for a terminal node @@ -148,23 +149,12 @@ NumConstNode::eval(const eval_context_type &eval_context) const throw (EvalExcep void NumConstNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const { - //CompileCode.write(reinterpret_cast<char *>(&FLDT), sizeof(FLDT)); - /*temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this)); - if (it != temporary_terms.end()) - { - CompileCode.write(&FLDT, sizeof(FLDT)); - idl= - CompileCode.write(reinterpret_cast<char *>(&idl),sizeof(idl)); - } - else - {*/ - CompileCode.write(&FLDC, sizeof(FLDC)); - double vard=atof(datatree.num_constants.get(id).c_str()); + CompileCode.write(&FLDC, sizeof(FLDC)); + double vard=atof(datatree.num_constants.get(id).c_str()); #ifdef DEBUGC - cout << "FLDC " << vard << "\n"; + cout << "FLDC " << vard << "\n"; #endif - CompileCode.write(reinterpret_cast<char *>(&vard),sizeof(vard)); - /*}*/ + CompileCode.write(reinterpret_cast<char *>(&vard),sizeof(vard)); } @@ -427,15 +417,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, double VariableNode::eval(const eval_context_type &eval_context) const throw (EvalException) { - // ModelTree::evaluateJacobian need to have the initval values applied to lead/lagged variables also - /*if (lag != 0) - throw EvalException();*/ - /*if(type==eModelLocalVariable) - cout << "eModelLocalVariable = " << symb_id << "\n";*/ eval_context_type::const_iterator it = eval_context.find(make_pair(symb_id, type)); if (it == eval_context.end()) { - //cout << "unknonw variable type = " << type << " simb_id = " << symb_id << "\n"; throw EvalException(); } @@ -445,15 +429,6 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep void VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const { - // If node is a temporary term - /*temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this)); - if (it != temporary_terms.end()) - { - CompileCode.write(&FLDT, sizeof(FLDT)); - int var=temporary_terms.count(const_cast<VariableNode *>(this))-1; - CompileCode.write(reinterpret_cast<char *>(&var), sizeof(var)); - return; - }*/ int i, lagl; #ifdef DEBUGC cout << "output_type=" << output_type << "\n"; @@ -501,6 +476,23 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou } } + +void +VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count, + temporary_terms_type &temporary_terms, + map<NodeID, pair<int, int> > &first_occurence, + int Curr_block, + Model_Block *ModelBlock, + int equation, + map_idx_type &map_idx) const +{ + if(type== eModelLocalVariable) + datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); +} + + + + void VariableNode::collectEndogenous(set<pair<int, int> > &result) const { @@ -711,9 +703,10 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, void UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const { NodeID this2 = const_cast<UnaryOpNode *>(this); @@ -721,8 +714,8 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, if (it == reference_count.end()) { reference_count[this2] = 1; - first_occurence[this2] = Curr_block; - arg->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + first_occurence[this2] = make_pair(Curr_block,equation); + arg->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); } else { @@ -730,7 +723,7 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) { temporary_terms.insert(this2); - ModelBlock->Block_List[first_occurence[this2]].Temporary_terms->insert(this2); + ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2); } } } @@ -1151,9 +1144,10 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, void BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const { NodeID this2 = const_cast<BinaryOpNode *>(this); @@ -1161,17 +1155,19 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, if (it == reference_count.end()) { reference_count[this2] = 1; - first_occurence[this2] = Curr_block; - arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); - arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + first_occurence[this2] = make_pair(Curr_block, equation); + arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); + arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); } else { reference_count[this2]++; if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) { + if(this2->idx==2280) + cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n"; temporary_terms.insert(this2); - ModelBlock->Block_List[first_occurence[this2]].Temporary_terms->insert(this2); + ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2); } } } @@ -1563,9 +1559,10 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, void TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const { NodeID this2 = const_cast<TrinaryOpNode *>(this); @@ -1573,18 +1570,20 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count, if (it == reference_count.end()) { reference_count[this2] = 1; - first_occurence[this2] = Curr_block; - arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); - arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); - arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, map_idx); + first_occurence[this2] = make_pair(Curr_block,equation); + arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); + arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); + arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx); } else { reference_count[this2]++; if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C) { + if(this2->idx==2280) + cout << "==>Curr_block= " << Curr_block << " equation= " << equation << " first_occurence[this2].first=" << first_occurence[this2].first << " first_occurence[this2].second=" << first_occurence[this2].second << "\n"; temporary_terms.insert(this2); - ModelBlock->Block_List[first_occurence[this2]].Temporary_terms->insert(this2); + ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2); } } } @@ -1743,9 +1742,10 @@ void UnknownFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output void UnknownFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const { cerr << "UnknownFunctionNode::computeTemporaryTerms: not implemented" << endl; diff --git a/MatlabFile.cc b/MatlabFile.cc new file mode 100644 index 0000000000000000000000000000000000000000..7fff2b498dd0565086d96769fd7b09b76ba18bc1 --- /dev/null +++ b/MatlabFile.cc @@ -0,0 +1,1076 @@ +/* + * Copyright (C) 2006-2008 Dynare Team + * + * This file is part of Dynare. + * + * Dynare 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. + * + * Dynare 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 Dynare. If not, see <http://www.gnu.org/licenses/>. + */ +/* +Usefull documentation: Matlab 7 Mat-File Format +----------------------------------------------- +revision: October 2008 PDF only Rereleased for Version 7.7 (Release 2008b) +available at: http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/matfile_format.pdf +*/ +#include "MatlabFile.hh" + + +MatlabFile::MatlabFile() +{ +} + +MatlabFile::~MatlabFile() +{ +} + +ArrayElem::ArrayElem() +{ +} + +ArrayElem::~ArrayElem() +{ +} + + +string +UTF8::ReadAlph(char* InBuff, int* pBuff, int Size) const +{ + char tmp_c[Size+1]; + string str; + memcpy(&tmp_c, InBuff+*pBuff, Size); + tmp_c[Size]=0; + str.assign(tmp_c); + if(Size<=4) + *pBuff += 4; + else if(Size % 8) + *pBuff += 8*ceil(double(Size) / 8); + else + *pBuff += Size; + return(str); +} + + +string +UTF16::ReadAlph(char* InBuff, int* pBuff, int Size) const +{ + string str(""); + for(int i=0;i<Size;i++) + str.append(&(InBuff[i*2+*pBuff])); + if(Size*2<=4) + *pBuff += 4; + else if((Size*2) % 8) + *pBuff += 8*ceil(double(Size*2) / 8); + else + *pBuff += Size*2; + return(str); +} + +string +UTF32::ReadAlph(char* InBuff, int* pBuff, int Size) const +{ + string str(""); + for(int i=0;i<Size;i++) + str.append(&(InBuff[i*4+*pBuff])); + if((Size*4) % 8) + *pBuff += 8*ceil(double(Size*4) / 8); + else + *pBuff += Size*4; + return(str); +} + +double +INT8::ReadNum(char* InBuff, int* pBuff) const +{ + char val; + val = InBuff[*pBuff]; + *pBuff += sizeof(val); + return(val); +} + + +double +INT16::ReadNum(char* InBuff, int* pBuff) const +{ + Int val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + +double +INT32::ReadNum(char* InBuff, int* pBuff) const +{ + LongInt val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + + +double +INT64::ReadNum(char* InBuff, int* pBuff) const +{ + LongLongInt val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + +double +Single::ReadNum(char* InBuff, int* pBuff) const +{ + float val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + +double +Double::ReadNum(char* InBuff, int* pBuff) const +{ + double val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + + + +Data_Header_t +SimpleElem::ReadDataHeader(char* InBuff, int* pBuff) const +{ + Data_Header_t data_header; + memcpy(&data_header.DataType, InBuff+*pBuff, sizeof(data_header.DataType)); + *pBuff += 2; + memcpy(&data_header.S_Number_of_Bytes, InBuff+*pBuff, sizeof(data_header.S_Number_of_Bytes)); + *pBuff += 2; + if(data_header.S_Number_of_Bytes!=0) + data_header.Number_of_Bytes=data_header.S_Number_of_Bytes; + else + { + memcpy(&data_header.Number_of_Bytes, InBuff+*pBuff, sizeof(data_header.Number_of_Bytes)); + *pBuff += 4; + } + return(data_header); +} + + +SimpleElem::SimpleElem() +{ + verbose=true; + array_elem=NULL; + Type=-1; +} + +SimpleElem::~SimpleElem() +{ +} + +returned_ReadData_t +SimpleElem::Get_Data_Class(Data_Header_t data_header) const +{ + returned_ReadData_t ret; + switch(data_header.DataType) + { + case miINT8: //8 bit, signed + case miUINT8: //8 bit, unsigned + ret.Type = Numerical; + ret.Simple = new INT8; + return(ret); + break; + case miINT16: //16-bit, signed + case miUINT16: //16-bit, unsigned + ret.Type = Numerical; + ret.Simple = new INT16; + return(ret); + break; + case miINT32: //32-bit, signed + case miUINT32: //32-bit, unsigned + ret.Type = Numerical; + ret.Simple = new INT32; + return(ret); + break; + case miINT64: //64-bit, signed + case miUINT64: //64-bit, unsigned + ret.Type = Numerical; + ret.Simple = new INT64; + return(ret); + break; + case miSINGLE: //IEEE� 754 single format + ret.Type = Numerical; + ret.Simple = new Single; + return(ret); + break; + case miDOUBLE: //IEEE 754 double format + ret.Type = Numerical; + ret.Simple = new Double; + return(ret); + break; + case miMATRIX: //MATLAB array + ret.Type = Matrix; + ret.Simple = NULL; + return(ret); + break; + case miCOMPRESSED: //Compressed Data + ret.Type = Compressed; + ret.Simple = NULL; + return(ret); + break; + case miUTF8: //Unicode UTF-8 Encoded Character Data + ret.Type = AlphaNumeric; + ret.Simple = new UTF8; + return(ret); + break; + case miUTF16: //Unicode UTF-16 Encoded Character Data + ret.Type = AlphaNumeric; + ret.Simple = new UTF16; + return(ret); + break; + case miUTF32: //Unicode UTF-32 Encoded Character Data + ret.Type = AlphaNumeric; + ret.Simple = new UTF32; + return(ret); + break; + default: + ret.Type = Unknown; + ret.Simple = NULL; + return(ret); + } +} + +void +SimpleElem::DataProceed(Data_Header data_header, char* InBuff, int* pBuff, FlagStructure_t flag) +{ + ArrayElem matrix; + returned_ReadData ret; + ret = Get_Data_Class(data_header); + Type = ret.Type; + double tmpv; + switch(ret.Type) + { + case Numerical: + if(data_header.Number_of_Bytes/ret.Simple->size()) + for(unsigned int i=0;i<data_header.Number_of_Bytes/ret.Simple->size();i++) + { + tmpv=ret.Simple->ReadNum(InBuff, pBuff); + VNumeric.push_back(tmpv); + } + //to align pBuff on a 4 Bytes + if(*pBuff % 4) + *pBuff += 4-(*pBuff % 4); + delete ret.Simple; + break; + case AlphaNumeric: + for(unsigned int i=0;i<data_header.Number_of_Bytes/ret.Simple->size();i++) + Vstr.push_back(ret.Simple->ReadAlph(InBuff, pBuff, data_header.Number_of_Bytes)); + //to align pBuff on a 4 Bytes + if(*pBuff % 4) + *pBuff += 4-(*pBuff % 4); + delete ret.Simple; + break; + case Matrix: + array_elem = matrix.ReadArray_class(InBuff, pBuff); + array_elem->ReadArray(InBuff, pBuff, flag); + break; + case Compressed: + cerr << "Error: Compressed data in Mat-file not implemnted yet!\n set option -v6 when saving to Mat-file\n"; + exit(EXIT_FAILURE); + case Unknown: + cerr << "Error: Mat-file format use incomptible format with Matlab 7 specification!\n set option -v6 when saving to Mat-file\n"; + exit(EXIT_FAILURE); + } +} + +returned_ReadData_t +SimpleElem::ReadData(char* InBuff, int* pBuff) const +{ + Data_Header_t data_header; + data_header = ReadDataHeader(InBuff, pBuff); + return(Get_Data_Class(data_header)); +} + +void +SimpleElem::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ + if(VNumeric.size()) + { + if(found) + collect_struct.variable_double_name[collect_struct.tmp_name]=VNumeric; + } + else if(Vstr.size()) + { + if(found) + collect_struct.variable_string_name[collect_struct.tmp_name]=Vstr; + } + else + { + if(array_elem) + { + array_elem->Collect(name, found, collect_struct); + } + } +} + + +void +SimpleElem::Print() const +{ + if(VNumeric.size()) + { + for(vector<double>::const_iterator it=VNumeric.begin(); it!=VNumeric.end(); it++) + cout << " " << *it; + } + else if(Vstr.size()) + { + for(vector<string>::const_iterator it=Vstr.begin(); it!=Vstr.end(); it++) + cout << " " << *it; + } + else + { + if(array_elem) + array_elem->Print(); + } +} + +void +SimpleElem::Delete() const +{ + if(array_elem) + { + array_elem->Delete(); + delete array_elem; + } +} + +LongInt +ArrayElem::ReadINT32(char* InBuff, int* pBuff) const +{ + LongInt val; + memcpy(&val, InBuff+*pBuff, sizeof(val)); + *pBuff += sizeof(val); + return(val); +} + + +Array_Flag_t +ArrayElem::ReadArrayFlag(char* InBuff, int* pBuff) /*const*/ +{ + Array_Flag_t array_flag; + memcpy(&array_flag, InBuff+*pBuff, sizeof(array_flag)); + *pBuff += sizeof(array_flag); + array_complex = (bool)(array_flag.flag & 16); + array_global = array_flag.flag & 32; + array_logical = array_flag.flag & 64; + type = (Array_Type)array_flag.classe; + if(type==Sparse_array) + array_nzmax = array_flag.nzmax; + return(array_flag); +} + +void +ArrayElem::ReadArrayDimension(char* InBuff, int* pBuff) /*const*/ +{ + Data_Header_t data_header; + data_header = ReadDataHeader(InBuff, pBuff); + number_of_dimensions = data_header.Number_of_Bytes/4; + for(int i=0; i<number_of_dimensions; i++) + { + double tmp_d=ReadINT32(InBuff, pBuff); + dimension.push_back(tmp_d); + } + if(number_of_dimensions % 2) + *pBuff += sizeof(dimension[0]); +} + +void +ArrayElem::ReadStructureNames(char* InBuff, int* pBuff) +{ + Data_Header_t data_header, data_header_2; + LongInt Field_name_length; + data_header=ReadDataHeader(InBuff, pBuff); + Field_name_length=ReadINT32(InBuff, pBuff); + data_header_2=ReadDataHeader(InBuff, pBuff); + Structure_number = data_header_2.Number_of_Bytes/Field_name_length; + char tmp_c[Field_name_length]; + for(int i=0; i<Structure_number;i++) + { + memcpy(tmp_c, InBuff+*pBuff, Field_name_length); + *pBuff += Field_name_length; + string variable_name(tmp_c); + Structure_Elem_name.push_back(variable_name); + } +} + +void +ArrayElem::ReadArrayName(char* InBuff, int* pBuff) /*const*/ +{ + Data_Header_t data_header; + data_header = ReadDataHeader(InBuff, pBuff); + char tmp_c[data_header.Number_of_Bytes+1]; + memcpy(&tmp_c, InBuff+*pBuff, data_header.Number_of_Bytes); + tmp_c[data_header.Number_of_Bytes]=0; + variable_name.assign(tmp_c); + if(data_header.Number_of_Bytes<=4) + *pBuff += 4; + else if(data_header.Number_of_Bytes % 8) + *pBuff += 8*ceil(double(data_header.Number_of_Bytes) / 8); + else + *pBuff += data_header.Number_of_Bytes; +} + + +ArrayElem* +ArrayElem::ReadArray_class(char* InBuff, int* pBuff) const +{ + Array_Flag_t array_flag; + ArrayElem array_elem; + array_flag=array_elem.ReadArrayFlag(InBuff, pBuff); + switch(array_flag.classe) + { + case Cell_array: + return(new CellArray); + break; + case Structure_: + return(new Structure); + break; + case Object_: + return(new Object); + break; + case Character_array: + return(new CharacterArray); + break; + case Sparse_array: + return(new SparseArray); + break; + case Double_precision_array: + return(new DoublePrecisionArray); + break; + case Single_precision_array: + return(new SinglePrecisionArray); + break; + case Signed_integer_8_bit: + return(new Bit8SignedInteger); + break; + case Unsigned_integer_8_bit: + return(new Bit8UnsignedInteger); + break; + case Signed_integer_16_bit: + return(new Bit16SignedInteger); + break; + case Unsigned_integer_16_bit: + return(new Bit16UnsignedInteger); + break; + case Signed_integer_32_bit: + return(new Bit32SignedInteger); + break; + case Unsigned_integer_32_bit: + return(new Bit32UnsignedInteger); + break; + default: + return(NULL); + } + return(NULL); +} + +void +ArrayElem::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +ArrayElem::Print() const +{ +} + +void +ArrayElem::Delete() const +{ +} + +void +CellArray::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +CellArray::Print() const +{ + //cout << "CellArray: "<< variable_name << "\n"; +} + +void +CellArray::Delete() const +{ + //cout << "CellArray: "<< variable_name << "\n"; +} + +void +CellArray::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + SimpleElem* simple; + ReadArrayDimension(InBuff, pBuff); + Cell_number = 1; + for(unsigned int i=0;i<dimension.size();i++) + Cell_number *= dimension[i]; + ReadArrayName(InBuff, pBuff); + flag.character=true; + for(int i=0;i<Cell_number;i++) + { + simple=new SimpleElem; + Data_Header_t data_header=simple->ReadDataHeader(InBuff, pBuff); + VCell.push_back(simple); + simple->DataProceed(data_header, InBuff, pBuff, flag); + } +} + + +void +Structure::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ + if(name==variable_name || found) + { + found = true; + vector<string>::const_iterator it2=Structure_Elem_name.begin(); + for(vector<PSimpleElem>::const_iterator it=VCell.begin(); it!=VCell.end(); it++) + { + collect_struct.tmp_name = *it2; + it2++; + (*it)->Collect(name, found, collect_struct); + } + } +} + +void +Structure::Print() const +{ + cout << "Structure: " << variable_name << "\n"; + vector<string>::const_iterator it2=Structure_Elem_name.begin(); + int i=0; + for(vector<PSimpleElem>::const_iterator it=VCell.begin(); it!=VCell.end(); it++) + { + cout << ++i << " -> " << *it2 << " :"; + it2++; + (*it)->Print(); + cout << "\n"; + } +} + +void +Structure::Delete() const +{ + vector<string>::const_iterator it2=Structure_Elem_name.begin(); + for(vector<PSimpleElem>::const_iterator it=VCell.begin(); it!=VCell.end(); it++) + { + (*it)->Delete(); + delete *it; + } +} + +void +Structure::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + SimpleElem* simple; + int i; + ReadArrayDimension(InBuff, pBuff); + ReadArrayName(InBuff, pBuff); + ReadStructureNames(InBuff, pBuff); + flag.no_name=true; + for(i=0;i<Structure_number;i++) + { + simple=new SimpleElem; + Data_Header_t data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); + data_header=simple->ReadDataHeader(InBuff, pBuff); + VCell.push_back(simple); + simple->DataProceed(data_header, InBuff, pBuff, flag); + } +} + +void +Object::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + + +void +Object::Print() const +{ +} + +void +Object::Delete() const +{ +} + + +void +Object::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + cerr << "Error: Object not implemented\n"; + exit(EXIT_FAILURE); +} + +void +CharacterArray::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ + if(name==variable_name || found) + { + found = true; + collect_struct.tmp_name=variable_name; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Collect(name, found, collect_struct); + } +} + +void +CharacterArray::Print() const +{ + cout << "CharacterArray: " << variable_name << "\n"; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Print(); + cout << "\n"; +} + +void +CharacterArray::Delete() const +{ + //cout << "CharacterArray: " << variable_name << "\n"; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Delete(); + delete *it; + //cout << "\n"; +} + + +void +CharacterArray::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + ReadArrayDimension(InBuff, pBuff); + Matrix_Elem_number = 1; + for(unsigned int i=0;i<dimension.size();i++) + { + Matrix_Elem_number *= dimension[i]; + } + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +SparseArray::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +SparseArray::Print() const +{ + cout << "Sparse Array: " << variable_name << "\n"; +} + +void +SparseArray::Delete() const +{ + //cout << "Sparse Array: " << variable_name << "\n"; +} + + +void +SparseArray::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + cerr << "Error: Sparse Array not implemented\n"; + exit(EXIT_FAILURE); +} + +void +DoublePrecisionArray::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ + if(name==variable_name || found) + { + found = true; + collect_struct.tmp_name=variable_name; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Collect(name, found, collect_struct); + } +} + + +void +DoublePrecisionArray::Print() const +{ + cout << "DoublePrecisionArray: " << variable_name << "\n"; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Print(); + cout << "\n"; +} + +void +DoublePrecisionArray::Delete() const +{ + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Delete(); + delete *it; +} + +void +DoublePrecisionArray::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + ReadArrayDimension(InBuff, pBuff); + Matrix_Elem_number = 1; + for(unsigned int i=0;i<dimension.size();i++) + { + Matrix_Elem_number *= dimension[i]; + } + if(!flag.no_name) + ReadArrayName(InBuff, pBuff); + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +SinglePrecisionArray::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ + if(name==variable_name || found) + { + found = true; + collect_struct.tmp_name=variable_name; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Collect(name, found, collect_struct); + } +} + + +void +SinglePrecisionArray::Print() const +{ + cout << "SinglePrecisionArray: " << variable_name << "\n"; + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Print(); + cout << "\n"; +} + +void +SinglePrecisionArray::Delete() const +{ + vector<PSimpleElem>::const_iterator it=VCell.begin(); + (*it)->Delete(); + delete *it; +} + +void +SinglePrecisionArray::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + ReadArrayDimension(InBuff, pBuff); + Matrix_Elem_number = 1; + for(unsigned int i=0;i<dimension.size();i++) + { + Matrix_Elem_number *= dimension[i]; + } + if(!flag.no_name) + ReadArrayName(InBuff, pBuff); + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit8SignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + + + +void +Bit8SignedInteger::Print() const +{ + //cout << "Bit8SignedInteger: \n"; +} + +void +Bit8SignedInteger::Delete() const +{ + //cout << "Bit8SignedInteger: \n"; +} + +void +Bit8SignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit8UnsignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + + +void +Bit8UnsignedInteger::Print() const +{ + //cout << "Bit8UnsignedInteger: \n"; +} + +void +Bit8UnsignedInteger::Delete() const +{ + //cout << "Bit8UnsignedInteger: \n"; +} + +void +Bit8UnsignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit16SignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +Bit16SignedInteger::Print() const +{ + //cout << "Bit16SignedInteger: \n"; +} + +void +Bit16SignedInteger::Delete() const +{ + //cout << "Bit16SignedInteger: \n"; +} + +void +Bit16SignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit16UnsignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +Bit16UnsignedInteger::Print() const +{ + //cout << "Bit16UnsignedInteger: \n"; +} + +void +Bit16UnsignedInteger::Delete() const +{ + //cout << "Bit16UnsignedInteger: \n"; +} + +void +Bit16UnsignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit32SignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +Bit32SignedInteger::Print() const +{ + //cout << "Bit32SignedInteger: \n"; +} + +void +Bit32SignedInteger::Delete() const +{ + //cout << "Bit32SignedInteger: \n"; +} + +void +Bit32SignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +void +Bit32UnsignedInteger::Collect(const string &name, bool found, CollectStruct &collect_struct) const +{ +} + +void +Bit32UnsignedInteger::Print() const +{ + //cout << "Bit32UnsignedInteger: \n"; +} + +void +Bit32UnsignedInteger::Delete() const +{ + //cout << "Bit32UnsignedInteger: \n"; +} + +void +Bit32UnsignedInteger::ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) +{ + Data_Header_t data_header; + SimpleElem* simple; + simple=new SimpleElem; + VCell.push_back(simple); + data_header=simple->ReadDataHeader(InBuff, pBuff); + simple->DataProceed(data_header, InBuff, pBuff, flag); +} + +Data_Header_t +MatlabFile::ReadDataHeader(ifstream &MatFile) +{ + Data_Header_t data_header; + MatFile.read(reinterpret_cast<char*>(&data_header.DataType),sizeof(data_header.DataType)); + MatFile.read(reinterpret_cast<char*>(&data_header.S_Number_of_Bytes),sizeof(data_header.S_Number_of_Bytes)); + if(data_header.S_Number_of_Bytes!=0) + data_header.Number_of_Bytes=data_header.S_Number_of_Bytes; + else + { + MatFile.read(reinterpret_cast<char*>(&data_header.Number_of_Bytes),sizeof(data_header.Number_of_Bytes)); + } + if(data_header.Number_of_Bytes<8) + data_header.Number_of_Bytes = 8; + return(data_header); +} + + + + + +void +MatlabFile::MatFileRead(string filename) +{ + ifstream MatFile; + Data_Header_t data_header; + SimpleElem *simpl; + int pBuff; + FlagStructure_t flag; + //ArrayElem elem; + MatFile.open(filename.c_str(),std::ios::in | std::ios::binary); + if (!MatFile.is_open()) + { + cerr << filename.c_str() << " Cannot be opened\n"; + exit(EXIT_FAILURE); + } + // Read the Header of the Mat-File + MatFile.read(reinterpret_cast<char*>(&header),sizeof(header)); + do + { + data_header=ReadDataHeader(MatFile); + char* InBuff; + InBuff = (char*)malloc(data_header.Number_of_Bytes+1); + MatFile.read(InBuff,data_header.Number_of_Bytes+1); + pBuff = 0; + simpl = new SimpleElem; + VSimpl.push_back(simpl); + flag.no_name=false; + flag.character=false; + simpl->DataProceed(data_header, InBuff, &pBuff, flag); + free(InBuff); + } + while(!MatFile.eof()); + MatFile.close(); +} + +bool +MatlabFile::Collect(const string &name, CollectStruct &collect_struct) const +{ + for(vector<PSimpleElem>::const_iterator it=VSimpl.begin(); it!=VSimpl.end(); it++) + (*it)->Collect(name, false, collect_struct); + return(!(collect_struct.variable_double_name.empty() and collect_struct.variable_string_name.empty())); +} + + +void +MatlabFile::MatFilePrint() +{ + for(vector<PSimpleElem>::iterator it=VSimpl.begin(); it!=VSimpl.end(); it++) + (*it)->Print(); +} + + +void +MatlabFile::Delete() +{ + for(vector<PSimpleElem>::iterator it=VSimpl.begin(); it!=VSimpl.end(); it++) + { + (*it)->Delete(); + delete *it; + } +} + +/* +int +main(int argc, char** argv) +{ + CollectStruct collect_struct; + MatlabFile matlab_file; + matlab_file.MatFileRead("gimf_steady.mat"); + //matlab_file.MatFileRead("essai.mat"); + matlab_file.MatFilePrint(); + bool tmp_b=false; + //string tmp_s("stored_values"); + tmp_b=matlab_file.Collect("stored_values", collect_struct); + if(tmp_b) + { + int i=0; + for(map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.begin();it2!=collect_struct.variable_double_name.end();it2++) + { + i++; + cout << i << " " << it2->first.c_str() << " : "; + for(vector<double>::iterator it=it2->second.begin();it!=it2->second.end();it++) + cout << *it; + cout << "\n"; + } + } +} + + +*/ diff --git a/ModFile.cc b/ModFile.cc index 554f80d76cfd822db1ef89f6a4b936e49e7e81cf..c2081049d8251bccfd8d21f69365c435ffaca68a 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -131,7 +131,6 @@ ModFile::evalAllExpressions() { if(global_eval_context.find(make_pair(j, eEndogenous))==global_eval_context.end()) { - //it2=mat_file.variable.find(symbol_table.getNameByID(eEndogenous, j)); map<string,vector<double> >::iterator it2=collect_struct.variable_double_name.find(symbol_table.getNameByID(eEndogenous, j)); if(it2!=collect_struct.variable_double_name.end()) { diff --git a/ModelTree.cc b/ModelTree.cc index 127d70578ef9c66f01b7b7917e8657d2127d6531..9c972fba235169d8ead26d32e09bdd2b68cbf328 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -277,7 +277,8 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) void ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) { - map<NodeID, int> reference_count, first_occurence; + map<NodeID, pair<int, int> > first_occurence; + map<NodeID, int> reference_count; int i, j, m, eq, var, lag; temporary_terms_type vect; ostringstream tmp_output; @@ -293,7 +294,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; - eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, map_idx); + eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx); } for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { @@ -303,7 +304,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i]; var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]; it=first_derivatives.find(make_pair(eq,variable_table.getID(eEndogenous, var,lag))); - it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, map_idx); + it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx); } } for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) @@ -314,7 +315,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; it=first_derivatives.find(make_pair(eq,variable_table.getID(eExogenous, var,lag))); - it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, map_idx); + it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx); } } //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr; @@ -328,7 +329,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock) eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; it=first_derivatives.find(make_pair(eq,variable_table.getID(eEndogenous, var,lag))); - it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, map_idx); + it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx); } } } @@ -470,7 +471,6 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & } output << " g2=0;g3=0;\n"; - temporary_terms_type tt2; if(ModelBlock->Block_List[j].Temporary_InUse->size()) { tmp_output.str(""); @@ -502,24 +502,24 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & sps = " "; else sps=""; - if (ModelBlock->Block_List[j].Temporary_terms->size()) - output << " " << sps << "% //Temporary variables" << endl; - i=0; - for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin(); - it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) - { - output << " " << sps; - (*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); - output << " = "; - (*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2); - // Insert current node into tt2 - tt2.insert(*it); - output << ";" << endl; - i++; - } // The equations for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { + temporary_terms_type tt2; + tt2.clear(); + if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size()) + output << " " << sps << "% //Temporary variables" << endl; + for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->begin(); + it != ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->end(); it++) + { + output << " " << sps; + (*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms); + output << " = "; + (*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2); + // Insert current node into tt2 + tt2.insert(*it); + output << ";" << endl; + } string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; lhs = eq_node->arg1; @@ -904,7 +904,6 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str << " % ////////////////////////////////////////////////////////////////////////" << endl; //The Temporary terms //output << global_output.str(); - temporary_terms_type tt2; if(ModelBlock->Block_List[j].Temporary_InUse->size()) { tmp_output.str(""); @@ -946,24 +945,24 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n"; } sps=""; - if (ModelBlock->Block_List[j].Temporary_terms->size()) - output << " " << sps << "% //Temporary variables" << endl; - i=0; - for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin(); - it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) - { - output << " " << sps; - (*it)->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); - output << " = "; - (*it)->writeOutput(output, oMatlabStaticModelSparse, tt2); - // Insert current node into tt2 - tt2.insert(*it); - output << ";" << endl; - i++; - } // The equations for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { + temporary_terms_type tt2; + tt2.clear(); + if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size()) + output << " " << sps << "% //Temporary variables" << endl; + for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->begin(); + it != ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->end(); it++) + { + output << " " << sps; + (*it)->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); + output << " = "; + (*it)->writeOutput(output, oMatlabStaticModelSparse, tt2); + // Insert current node into tt2 + tt2.insert(*it); + output << ";" << endl; + } string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ; output << sps << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl; @@ -1202,38 +1201,41 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl } else lhs_rhs_done=false; - //The Temporary terms - temporary_terms_type tt2; - i=0; - for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin(); - it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) + // The equations + for (i = 0;i < ModelBlock->Block_List[j].Size;i++) { - (*it)->compile(code_file,false, output_type, tt2, map_idx); - code_file.write(&FSTPT, sizeof(FSTPT)); - map_idx_type::const_iterator ii=map_idx.find((*it)->idx); - v=(int)ii->second; - code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); - // Insert current node into tt2 - tt2.insert(*it); + //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); + //The Temporary terms + temporary_terms_type tt2; #ifdef DEBUGC - cout << "FSTPT " << v << "\n"; - code_file.write(&FOK, sizeof(FOK)); - code_file.write(reinterpret_cast<char *>(&i), sizeof(i)); + k=0; #endif - i++; - } - for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin(); - it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) - { - map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->begin(); + it != ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->end(); it++) + { + (*it)->compile(code_file,false, output_type, tt2, map_idx); + code_file.write(&FSTPT, sizeof(FSTPT)); + map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + v=(int)ii->second; + code_file.write(reinterpret_cast<char *>(&v), sizeof(v)); + // Insert current node into tt2 + tt2.insert(*it); #ifdef DEBUGC - cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n"; + cout << "FSTPT " << v << "\n"; + code_file.write(&FOK, sizeof(FOK)); + code_file.write(reinterpret_cast<char *>(&k), sizeof(k)); + ki++; +#endif + + } +#ifdef DEBUGC + for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin(); + it != ModelBlock->Block_List[j].Temporary_terms->end(); it++) + { + map_idx_type::const_iterator ii=map_idx.find((*it)->idx); + cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n"; + } #endif - } - // The equations - for (i = 0;i < ModelBlock->Block_List[j].Size;i++) - { - //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0); if (!lhs_rhs_done) { eq_node = equations[ModelBlock->Block_List[j].Equation[i]]; diff --git a/include/ExprNode.hh b/include/ExprNode.hh index 06d9228e019d7b00174f82e7a04709a8bf594374..9fb52b9df2f1b4aabb57cf0fcf88c0fa774a1a7c 100644 --- a/include/ExprNode.hh +++ b/include/ExprNode.hh @@ -148,9 +148,10 @@ public: virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const = 0; virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const; class EvalException @@ -204,6 +205,13 @@ public: virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const; virtual void collectEndogenous(set<pair<int, int> > &result) const; virtual void collectExogenous(set<pair<int, int> > &result) const; + virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, + temporary_terms_type &temporary_terms, + map<NodeID, pair<int, int> > &first_occurence, + int Curr_block, + Model_Block *ModelBlock, + int equation, + map_idx_type &map_idx) const; virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const; virtual double eval(const eval_context_type &eval_context) const throw (EvalException); virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const; @@ -225,9 +233,10 @@ public: virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const; virtual void collectEndogenous(set<pair<int, int> > &result) const; virtual void collectExogenous(set<pair<int, int> > &result) const; @@ -254,9 +263,10 @@ public: virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const; virtual void collectEndogenous(set<pair<int, int> > &result) const; virtual void collectExogenous(set<pair<int, int> > &result) const; @@ -290,9 +300,10 @@ public: virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const; virtual void collectEndogenous(set<pair<int, int> > &result) const; virtual void collectExogenous(set<pair<int, int> > &result) const; @@ -317,9 +328,10 @@ public: virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, - map<NodeID, int> &first_occurence, + map<NodeID, pair<int, int> > &first_occurence, int Curr_block, Model_Block *ModelBlock, + int equation, map_idx_type &map_idx) const; virtual void collectEndogenous(set<pair<int, int> > &result) const; virtual void collectExogenous(set<pair<int, int> > &result) const; @@ -349,7 +361,8 @@ struct Block bool is_linear; int *Equation, *Own_Derivative; int *Variable, *Other_Endogenous, *Exogenous; - temporary_terms_type *Temporary_terms; + temporary_terms_type **Temporary_Terms_in_Equation; + //temporary_terms_type *Temporary_terms; temporary_terms_inuse_type *Temporary_InUse; IM_compact *IM_lead_lag; int Code_Start, Code_Length; diff --git a/include/MatlabFile.hh b/include/MatlabFile.hh new file mode 100644 index 0000000000000000000000000000000000000000..807b1743959776441666f896638ed45c63f9c673 --- /dev/null +++ b/include/MatlabFile.hh @@ -0,0 +1,426 @@ +/* + * Copyright (C) 2006-2008 Dynare Team + * + * This file is part of Dynare. + * + * Dynare 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. + * + * Dynare 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 Dynare. If not, see <http://www.gnu.org/licenses/>. + */ +/* +Usefull documentation: Matlab 7 Mat-File Format +----------------------------------------------- +revision: October 2008 PDF only Rereleased for Version 7.7 (Release 2008b) +available at: http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/matfile_format.pdf +*/ + +#ifndef _MAT_FILE_HH +#define _MAT_FILE_HH + +#include <cstdio> +#include <cstdlib> +#include <string> +#include <fstream> +#include <cstring> +#include <iostream> +#include <map> +#include <vector> +#include <math.h> +//! zlib needed to uncompress the mat-file. It is available with GCC 4.3.2 but it needs a dll !! +//! => to avoid compress MatFile, save is used with option '-v6' in save_params_and_strady_state.m +//#include "zlib.h" +using namespace std; + + +enum Data_Type +{ + miINT8 = 1, //8 bit, signed + miUINT8 = 2, //8 bit, unsigned + miINT16 = 3, //16-bit, signed + miUINT16 = 4, //16-bit, unsigned + miINT32 = 5, //32-bit, signed + miUINT32 = 6, //32-bit, unsigned + miSINGLE = 7, //IEEE� 754 single format + miDOUBLE = 9, //IEEE 754 double format + miINT64 = 12, //64-bit, signed + miUINT64 = 13, //64-bit, unsigned + miMATRIX = 14, //MATLAB array + miCOMPRESSED = 15, //Compressed Data + miUTF8 = 16, //Unicode UTF-8 Encoded Character Data + miUTF16 = 17, //Unicode UTF-16 Encoded Character Data + miUTF32 = 18 //Unicode UTF-32 Encoded Character Data +}; + +enum Array_Type +{ + Cell_array = 1, + Structure_ = 2, + Object_ = 3, + Character_array = 4, + Sparse_array = 5, + Double_precision_array = 6, + Single_precision_array = 7, + Signed_integer_8_bit = 8, + Unsigned_integer_8_bit = 9, + Signed_integer_16_bit = 10, + Unsigned_integer_16_bit = 11, + Signed_integer_32_bit = 12, + Unsigned_integer_32_bit = 13 +}; + + +enum Function_Returned_Type +{ + Numerical =1, + AlphaNumeric =2, + Matrix =3, + Compressed =4, + Unknown =5 +}; + +class ArrayElem; +class SimpleElem; + +typedef long long LongLongInt; +typedef long int LongInt; +typedef short int Int; +//typedef char ShortInt; +typedef unsigned long int uLongInt ; +typedef unsigned short int uShortInt ; +typedef short int ShortInt ; +typedef class SimpleElem *PSimpleElem; + +//!Header of MatFile +typedef struct Header +{ + char Theader[124]; + short int Version; + char Edian_Indicator[2]; +} +Header_t; + +typedef struct Data_Header +{ + ShortInt S_Number_of_Bytes; + ShortInt DataType; + uLongInt Number_of_Bytes; +} +Data_Header_t; + +typedef struct Array_Flag +{ + Data_Header_t tag; + unsigned char classe; + unsigned char flag; + char undef1[2]; + uLongInt nzmax; +} +Array_Flag_t; + + +typedef struct returned_ReadData +{ + SimpleElem* Simple; + Function_Returned_Type Type; +} returned_ReadData_t; + + +typedef struct FlagStructure +{ + bool no_name; + bool character; +} FlagStructure_t; + +typedef struct CollectStruct +{ + /*vector<string> variable_name; + vector< vector<double> > variable_double; + vector< vector<string> > variable_string;*/ + string tmp_name; + map<string,vector<string> > variable_string_name; + map<string,vector<double> > variable_double_name; +} +CollectStruct; + +typedef vector<int> Array_Dimensions_t; + + + + + +//! Base class for simple elements in Mat-File +class SimpleElem +{ +public: + bool verbose; + vector<double> VNumeric; + vector<string> Vstr; + ArrayElem *array_elem; + int Type; + SimpleElem(); + virtual ~SimpleElem(); + virtual double ReadNum(char* InBuff, int* pBuff) const{return(0);}; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const{return(NULL);}; + virtual Data_Header_t ReadDataHeader(char* InBuff, int* pBuff) const; + virtual int size() const{cout << "oups\n";return(0);}; + void Print() const; + void Delete() const; + void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + returned_ReadData_t ReadData(char* InBuff, int* pBuff) const; + returned_ReadData_t Get_Data_Class(Data_Header_t data_header) const; + void DataProceed(Data_Header data_header, char* InBuff, int* pBuff, FlagStructure flag); +}; + + +class UTF8 : public SimpleElem +{ +public: + virtual int size() const {return(1);}; + virtual double ReadNum(char* InBuff, int* pBuff) const {return(NULL);}; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const; +}; + +class UTF16 : public SimpleElem +{ +public: + virtual int size() const {return(2);}; + virtual double ReadNum(char* InBuff, int* pBuff) const {return(NULL);}; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const; +}; + +class UTF32 : public SimpleElem +{ +public: + virtual int size() const {return(4);}; + virtual double ReadNum(char* InBuff, int* pBuff) const {return(NULL);}; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const; +}; + +class INT8 : public SimpleElem +{ +public: + virtual int size() const {return(1);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + +class INT16 : public SimpleElem +{ +public: + virtual int size() const {return(2);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + + +class INT32 : public SimpleElem +{ +public: + virtual int size() const {return(4);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + + +class INT64 : public SimpleElem +{ +public: + virtual int size() const {return(8);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + + +class Single : public SimpleElem +{ +public: + virtual int size() const {return(4);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + + +class Double : public SimpleElem +{ +public: + virtual int size() const {return(8);}; + virtual double ReadNum(char* InBuff, int* pBuff) const; + virtual string ReadAlph(char* InBuff, int* pBuff, int Size) const {return(NULL);}; +}; + +//! Base class for Array Element in Mat-File +class ArrayElem : public SimpleElem +{ +private: + +protected: + +public: + int Cell_number, Structure_number, Matrix_Elem_number; + Array_Type Type; + vector<PSimpleElem> VCell; + vector<string> Structure_Elem_name; + bool array_complex, array_global, array_logical; + string variable_name; + int array_nzmax; + int number_of_dimensions; + vector<int> dimension; + vector<double> Double_value; + vector<string> String_value; + + Array_Type type; + ArrayElem(); + virtual LongInt ReadINT32(char* InBuff, int* pBuff) const; + virtual Array_Flag_t ReadArrayFlag(char* InBuff, int* pBuff) /*const*/; + virtual void ReadArrayDimension(char* InBuff, int* pBuff) /*const*/; + virtual void ReadArrayName(char* InBuff, int* pBuff) /*const*/; + virtual void ReadStructureNames(char* InBuff, int* pBuff); + virtual ArrayElem* ReadArray_class(char* InBuff, int* pBuff) const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag) /*const*/ {cout << "oups..\n";}; + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual ~ArrayElem(); +}; + +class CellArray : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Structure : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Object : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class CharacterArray : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class SparseArray : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class DoublePrecisionArray : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class SinglePrecisionArray : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit8SignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit8UnsignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit16SignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit16UnsignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit32SignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + +class Bit32UnsignedInteger : public ArrayElem +{ + public: + virtual void Collect(const string &name, bool found, CollectStruct &collect_struct) const; + virtual void Print() const; + virtual void Delete() const; + virtual void ReadArray(char* InBuff, int* pBuff, FlagStructure_t flag); +}; + + +class MatlabFile +{ + public: + Header_t header; + vector<PSimpleElem> VSimpl; + MatlabFile(); + ~MatlabFile(); + void MatFileRead(string filename); + void MatFilePrint(); + void Delete(); + bool Collect(const string &name, CollectStruct &collect_struct) const; + Data_Header_t ReadDataHeader(ifstream &MatFile); +}; +#endif