function.cpp 30.3 KB
Newer Older
1 2 3
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

   Copyright (C) 2015 CNRS
4
   Copyright (C) 2013, 2014, 2016 Inria
5 6 7 8 9 10 11

   This file is part of ALTA.

   This Source Code Form is subject to the terms of the Mozilla Public
   License, v. 2.0.  If a copy of the MPL was not distributed with this
   file, You can obtain one at http://mozilla.org/MPL/2.0/.  */

12 13 14
#include "function.h"		

#include "common.h"
15
#include "params.h"
16 17
#include "plugins_manager.h"

18 19
#include <cassert>

20 21
using namespace alta;

22 23
/*--- Functions implementation ----*/

24
void function::bootstrap(const ptr<data>, const arguments& args)
25
{
26 27 28 29 30 31 32 33 34
  #ifdef BOOTSTRAP_DEBUG
  std::cout << __FILE__ << " " << __LINE__ << " args = " << args << std::endl;
  #endif

  // If the bootstrap option contains a filename, load it
  if(args.is_defined("bootstrap") && !args["bootstrap"].empty())
  {
    std::ifstream in(args["bootstrap"].c_str());
    if(in.is_open())
35
    {
36 37 38 39 40 41 42 43 44
      if(!load(in))
			{
				std::cerr << "<<ERROR>> cannot bootstrap from file \"" << args["bootstrap"] << "\"" << std::endl;
			}
    }
    else // No file for the boostrap arguments
    {
      std::cerr << " Bootstrapping from command line is not implemented" << std::endl;
			NOT_IMPLEMENTED();
45
    }
46 47
  }//end of if bootstrap is defined and not empty
    
48 49
}

50 51
void function::save(const std::string& filename, const arguments& args) const
{
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
52 53 54 55 56
	bool const is_alta     = !args.is_defined("export") || args["export"] == "alta";
	bool const is_cpp      = args["export"] == "C++";
	bool const is_explorer = args["export"] == "explorer";
	bool const is_shader   = args["export"] == "shader" || is_explorer;
	bool const is_matlab   = args["export"] == "matlab";
57

58
	// Open the file
59
	std::ofstream file(filename.c_str());
60 61 62 63 64
	if(!file.is_open())
	{
		std::cerr << "<<ERROR>> unable to open output file for writing" << std::endl;
	}

65 66 67 68 69 70 71 72 73 74 75 76
	// If the export is the alta format, use the maximum precision formatting
	if(is_alta)
	{
		file.precision(10);
	}

	if(is_explorer)
	{
		file << "analytic" << std::endl;
		file << std::endl;
		file << "::begin shader" << std::endl;
	}
77

78
	// Save common header
79
	save_header(file, args);
80 81

	// Save function definition
82
	save_body(file, args);
83

84 85 86 87
	if(is_cpp)
	{
		file << "vec brdf(const vec& in, const vec& file)" << std::endl;
		file << "{" << std::endl;
88
		file << "\tvec res(" << parametrization().dimY() << ");" << std::endl;
89 90 91 92 93 94
		file << "\t";
	}
	else if(is_matlab)
	{
		file << "function res = brdf(in, file)" << std::endl;
		file << "{" << std::endl;
95
		file << "\tres = zeros(" << parametrization().dimY() << ");" << std::endl;
96 97 98 99 100 101 102 103
		file << "\t";
	}
	else if(is_shader)
	{
		file << "vec3 BRDF(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y)" << std::endl;
		file << "{" << std::endl;
		file << "\tvec3 res = ";
	}
104 105

	// Save fit data
106 107 108 109 110 111 112 113 114 115 116 117 118 119
	save_call(file, args);
	if(is_cpp || is_shader)
	{
		file << ";" << std::endl;
		file << "\treturn res;" << std::endl;
		file << "}" << std::endl;
	}
	else if(is_matlab)
	{
		file << ";" << std::endl;
		file << "\treturn res;" << std::endl;
		file << "endfunction" << std::endl;
	}
	file << std::endl;
120

121 122 123 124
	if(is_explorer)
	{
		file << "::end shader" << std::endl;
	}
125 126 127 128 129 130 131 132 133
}
		
//! \brief save the header of the output function file. The header should
//! store general information about the fit such as the command line used
//! the dimension of the fit. L2 and L_inf distance could be added here.
void function::save_header(std::ostream& out, const arguments& args) const
{
	if(!args.is_defined("export"))
	{
134
		out << "#ALTA FUNC HEADER" << std::endl;
135
		out << "#CMD " << args.get_cmd() << std::endl;
136 137
		out << "#DIM " << parametrization().dimX() << " " << parametrization().dimY() << std::endl;
		out << "#PARAM_IN  " << params::get_name(parametrization().input_parametrization()) << std::endl;
138
		//out << "#PARAM_OUT " << params::get_name(output_parametrization()) << std::endl;*
139 140 141 142
		if(args.is_defined("export-append")) 
		{
			out << args["export-append"] << std::endl;
		}
143
		out << "#ALTA HEADER END" << std::endl;
144 145 146 147 148 149 150
		out << std::endl;
	}
}

//! \brief save function specific data. This has no use for ALTA export
//! but allows to factorize the code in the C++ or matlab export by
//! defining function calls that are common to all the plugins.
151
void function::save_body(std::ostream&, const arguments&) const
152 153
{

154

155 156 157 158 159
}

//! \brief save object specific information. For an ALTA export the
//! coefficients will be exported. For a C++ or matlab export, the call
//! to the associated function will be done.
160
void function::save_call(std::ostream&, const arguments&) const
161 162 163 164
{

}

165
//! \brief L2 norm to data.
166
double function::L2_distance(const ptr<data>& d) const
167
{
168
	double l2_dist = 0.0;
169 170 171
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
172 173
    vec x(parametrization().dimX());
    vec y(parametrization().dimY());
174

175
		if(parametrization().input_parametrization() == params::UNKNOWN_INPUT)
176
		{
177
			memcpy(&x[0], &dat[0], parametrization().dimX()*sizeof(double));
178 179 180
		}
		else
		{
181 182
      params::convert(&dat[0],
                      d->parametrization().input_parametrization(),
183
                      parametrization().input_parametrization(), &x[0]);
184
		}
185 186
    memcpy(&y[0], &dat[d->parametrization().dimX()],
           parametrization().dimY()*sizeof(double));
187

PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
188
    l2_dist += std::pow(norm(y-value(x)), 2);
189
	}
190 191

  assert(d->size() > 0.);
192
	l2_dist = std::sqrt(l2_dist / d->size());
193

194
	return l2_dist;
195 196 197
}

//! \brief Linf norm to data.
198
double function::Linf_distance(const ptr<data>& d) const
199
{
200 201
	vec mean = vec::Zero(parametrization().dimY());
	vec var  = vec::Zero(parametrization().dimY());
202

203
#ifdef DEBUG
204
	std::cout << "<<DEBUG>> input param here = " << params::get_name(parametrization().input_parametrization()) << std::endl;
205
#endif
206

207
	double linf_dist = 0.0;
208 209 210
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
211
		vec x(parametrization().dimX()), y(parametrization().dimY());
212

PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
213 214
    // Convert the position of the data sample to the parametrization
    // of the function.
215
    if(parametrization().input_parametrization() == params::UNKNOWN_INPUT)
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
216
    {
217
        memcpy(&x[0], &dat[0], parametrization().dimX()*sizeof(double));
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
218 219 220
    }
    else
    {
221 222
        params::convert(&dat[0],
                        d->parametrization().input_parametrization(),
223
                        parametrization().input_parametrization(), &x[0]);
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
224
    }
225 226 227

		// Copy the value part of the data vector in a vector to perform vector
		// operations on it (used in the computation of the mean).
228 229
    memcpy(&y[0], &dat[d->parametrization().dimX()],
           d->parametrization().dimY() * sizeof(double));
230 231 232

		// Take the componentwise-max of the two vectors.
		const vec v = value(x);
233
    for(int j=0; j<d->parametrization().dimY(); ++j)
234
		{
235
			linf_dist = std::max<double>(linf_dist, std::abs(y[j]-v[j]));
236
		}
237

238
		// Compute the mean
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
239
		mean += (y-v) / static_cast<double>(d->size());
240 241
	}
	
242
	// Compute the standard deviation with respect to the mean error
243 244 245
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
246
    vec x(parametrization().dimX()), y(d->parametrization().dimY()), val(parametrization().dimY());
247 248 249
		
        // Convert the position of the data sample to the parametrization
        // of the function.
250
        if(parametrization().input_parametrization() == params::UNKNOWN_INPUT)
251
        {
252
            memcpy(&x[0], &dat[0], parametrization().dimX()*sizeof(double));
253 254 255
        }
        else
        {
256 257
            params::convert(&dat[0],
                            d->parametrization().input_parametrization(),
258
                            parametrization().input_parametrization(),
259
                            &x[0]);
260
        }
261

262 263
		// Copy the value part of the data vector in a vector to perform vector
		// operations on it (used in the computation of the mean).
264
    memcpy(&y[0], &dat[d->parametrization().dimX()], parametrization().dimY()*sizeof(double));
265

266
		val = value(x);
267
    for(int j=0; j<d->parametrization().dimY(); ++j)
268
		{ 
269
      y[j] = dat[d->parametrization().dimX()+j];
270 271
			var[j] += pow(mean[j] - (val[j]-y[j]), 2) / double(d->size());
		}
272 273
	}

274 275
	std::cout << "<<INFO>> Mean = " << mean << ", Var = " << var << std::endl;

276
	return linf_dist;
277
}
278

279

280
/*--- Non-linear functions implementation ----*/
281 282 283 284 285

nonlinear_function::nonlinear_function()
{
}

286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
bool nonlinear_function::load(std::istream& in)
{
	// Parse line until the next comment
	while(in.peek() != '#')
	{
		char line[256];
		in.getline(line, 256);

		// If we cross the end of the file, or the badbit is
		// set, the file cannot be loaded
		if(!in.good())
			return false;
	}

	// Checking for the comment line #FUNC nonlinear_function_phong
	std::string token;
	in >> token;
	if(token.compare("#FUNC") != 0) 
	{ 
		std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl; 
#ifdef DEBUG
		std::cout << "<<DEBUG>> got: \"" << token << "\"" << std::endl;
#endif
		return false;
	}

	in >> token;
	if(token.compare("nonlinear_function") != 0)
	{
		std::cerr << "<<ERROR>> parsing the stream. A function name is defined." << std::endl;
		std::cerr << "<<ERROR>> did you forget to specify the plugin used to export?" << std::endl;
		return false;
	}

	int nb_params = nbParameters();
	vec p(nb_params);
	for(int i=0; i<nb_params; ++i)
	{
		in >> token >> p[i];
	}

	setParameters(p);
	return true;
}

void nonlinear_function::save_call(std::ostream& out, const arguments& args) const
{
	if(!args.is_defined("export"))
	{
		// Dump a #FUNC nonlinear
		out << "#FUNC nonlinear_function" << std::endl;

		// Dump the parameters in order
		vec p = parameters();
		for(int i=0; i<p.size(); ++i)
		{
			out << "param_" << i+1 << "\t" << p[i] << std::endl;
		}
		out << std::endl;
	}

	function::save_call(out, args);
}

350
void nonlinear_function::bootstrap(const ptr<data> d, const arguments& args)
351
{
352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
  #ifdef BOOTSTRAP_DEBUG
  std::cout << __FILE__ << " " << __LINE__ << " args = " << args << std::endl;
  #endif

  if(args.is_vec("bootstrap"))
  {
    vec p = args.get_vec("bootstrap", nbParameters());

    #ifdef BOOTSTRAP_DEBUG
    std::cout << __FILE__ << " " << __LINE__ << " args = " << args << std::endl;
    std::cout << "BOOTSTRAPPING VALUE = " << p << std::endl
              << " with nbParameters = " << nbParameters() << std::endl;
    #endif

    setParameters(p);
  }
  else
  {
    #ifdef BOOTSTRAP_DEBUG
    std::cout << __FILE__ << " " << __LINE__ << " args = " << args << std::endl;
    #endif

    function::bootstrap(d, args);
  }
376
}
377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396
		
vec nonlinear_function::getParametersMax() const
{
	vec M(nbParameters());
	for(int i=0; i<nbParameters(); ++i) 
	{
		M[i] = std::numeric_limits<double>::max();
	}
	return M;
}

vec nonlinear_function::getParametersMin() const
{
	vec m(nbParameters());
	for(int i=0; i<nbParameters(); ++i) 
	{
		m[i] = -std::numeric_limits<double>::max();
	}
	return m;
}
397

398 399 400


/*--- Compound functions implementation ----*/
401 402 403
compound_function::~compound_function()
{
}
404 405 406 407 408 409 410

vec compound_function::operator()(const vec& x) const
{
	return value(x);
}
vec compound_function::value(const vec& x) const
{
411
	vec res = vec::Zero(parametrization().dimY());
412 413
	for(unsigned int i=0; i<fs.size(); ++i)
	{
414 415 416
		vec temp_x(fs[i]->parametrization().dimX());
    params::convert(&x[0], parametrization().input_parametrization(),
                    fs[i]->parametrization().input_parametrization(), &temp_x[0]);
417 418 419 420 421 422 423 424
		res = res + fs[i]->value(temp_x);
	}
	return res;
}

vec compound_function::parametersJacobian(const vec& x) const
{
	int nb_params = nbParameters();
425
	vec jac = vec::Zero(nb_params*parametrization().dimY());
426 427 428 429 430 431

	int start_i = 0;

	// Export the sub-Jacobian for each function
	for(unsigned int f=0; f<fs.size(); ++f)
	{
432
		const ptr<nonlinear_function>& func = fs[f];
433 434 435 436
		int nb_f_params = func->nbParameters(); 

		// Only export Jacobian if there are non-linear parameters
		if(nb_f_params > 0 && !is_fixed[f])
437
		{
438

439 440 441
			vec temp_x(func->parametrization().dimX());
      params::convert(&x[0], parametrization().input_parametrization(),
                      func->parametrization().input_parametrization(), &temp_x[0]);
442 443 444 445
			vec func_jac = func->parametersJacobian(temp_x);

			for(int i=0; i<nb_f_params; ++i)
			{
446
				for(int y=0; y<parametrization().dimY(); ++y)
447 448 449 450 451 452
				{
					jac[y*nb_params + (i+start_i)] = func_jac[y*nb_f_params + i];
				}
			}

			start_i += nb_f_params;
453 454
		}
	}
455 456

	return jac;
457 458
}

459 460 461
// Return the parameters of the composition of FUNCTIONS.
static const parameters
compound_parameters(const std::vector<ptr<nonlinear_function> >& functions)
462
{
463
    assert(functions.size() > 0);
464

465
    parameters result = functions[0]->parametrization();
466

467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488
    for (auto f: functions)
    {
        auto input = f->parametrization().input_parametrization();
        auto output = f->parametrization().output_parametrization();

        if (result.input_parametrization() == params::UNKNOWN_INPUT)
            result = alta::parameters(result.dimX(), result.dimY(),
                                      input,
                                      result.output_parametrization());
        else if (result.input_parametrization() != input)
            // Parametrization mismatch: fall back to Cartesian.
            result = result.set_input(6, params::CARTESIAN);

        if (result.output_parametrization() == params::UNKNOWN_OUTPUT)
            result = alta::parameters(result.dimX(), result.dimY(),
                                      result.input_parametrization(),
                                      output);
        else if (result.output_parametrization() != output)
            abort();
    }

    return result;
489 490
}

491 492 493 494 495 496 497 498 499 500 501 502 503 504 505
compound_function::compound_function(const std::vector<ptr<nonlinear_function> >& functions,
                                     const std::vector<arguments> args)
    : nonlinear_function(compound_parameters(functions)),
      fs(functions), fs_args(args)
{
    assert(functions.size() == args.size());

    is_fixed.resize(functions.size());
    for (unsigned int x = 0; x < functions.size(); x++)
    {
        is_fixed[x] = args[x].is_defined("fixed");
    }
}


506 507 508 509 510
nonlinear_function* compound_function::operator[](int i) const
{
#ifdef DEBUG
	assert(i >= 0 && i < fs.size());
#endif
511
	return fs[i].get();
512
}
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533
		
unsigned int compound_function::size() const
{
	return fs.size();
}

bool compound_function::load(std::istream& in)
{
	int nb_good = 0; // Number of correctly openned functions
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		std::streampos pos = in.tellg();
		if(!fs[i]->load(in))
		{
			in.seekg(pos);
		}
		else
		{
			nb_good++;
		}
	}
534
	std::cout << "<<DEBUG>> number of correct function loaded in compound: " << nb_good  << " / " << fs.size() << std::endl;
535 536 537
	return nb_good > 0;
}

538
void compound_function::bootstrap(const ::ptr<data> d, const arguments& args)
539
{
540 541 542 543 544
	bool const global_bootstrap = args.is_defined("bootstrap");

	#ifdef BOOTSTRAP_DEBUG
	std::cout << __FILE__ << " " << __LINE__ << " args = " << args << std::endl;
	#endif
545

546 547 548 549 550 551 552
	// Check if the bootstrap is global
	if(global_bootstrap)
	{
		if(args.is_vec("bootstrap"))
		{
			vec p = args.get_vec("bootstrap", nbParameters());
			setParameters(p);
553 554 555 556 557 558 559 560
		
			for(unsigned int i=0; i<fs.size(); ++i)
			{
				if(fs[i]->nbParameters() == 0)
				{
					fs[i]->bootstrap(d, fs_args[i]);
				}
			}
561 562
      std::cout << "<<DEBUG>> Number of parameters for compound function: " << nbParameters() << std::endl;
      std::cout << "<<INFO>> Will use " << p << " as a bootstrap vector for the compound non-linear function" << std::endl;
563 564 565 566 567 568 569 570 571 572 573 574 575 576
		}
		else
		{
			std::ifstream file;
			file.open(args["bootstrap"].c_str()) ;
			if(file.is_open())
			{
				// Skip the header
				std::string line ;
				do
				{
					std::getline(file, line) ;
				}
				while(line != "#ALTA HEADER END");
577

578 579 580 581 582 583 584
				// Load the file
				// For each individual element, the input file is parsed
				// in order. if the function cannot parse it, the file is
				// restored to the previous state and the loading continues
				for(unsigned int i=0; i<fs.size(); ++i)
				{
					std::streampos pos = file.tellg();
585
					
586 587 588
					ptr<product_function> p = dynamic_pointer_cast<product_function>(fs[i]);
					if(p) {

589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
						nonlinear_function* f1 = p->first();
						nonlinear_function* f2 = p->second();

						pos = file.tellg();
						if(!f1->load(file))
						{
							file.seekg(pos);

							// Bootstrap the function as if it was not loaded
							f1->bootstrap(d, args);
							std::cout << "<<DEBUG>> Unable to load first function of product, regular bootstraping" << std::endl;
						}

						// If the second function cannot be loaded, put the input stream
						// in the previous state and bootstrap normaly this function.
						pos = file.tellg();
						if(!f2->load(file))
						{
							file.seekg(pos);

							// Bootstrap the function as if it was not loaded
							f2->bootstrap(d, args);
							std::cout << "<<DEBUG>> Unable to load second function of product, regular bootstraping" << std::endl;
						}
613
					}
614 615 616 617 618 619 620
					else
					{	
						// If the function cannot be loaded, put the input stream
						// in the previous state and bootstrap normaly this function.
						if(!fs[i]->load(file))
						{
							file.seekg(pos);
621

622 623 624 625
							// Bootstrap the function as if it was not loaded
							fs[i]->bootstrap(d, fs_args[i]);
							std::cout << "<<DEBUG>> Unable to load one function of compound, regular bootstraping" << std::endl;
						}
626 627 628
					}//end of else
				}//end of for-loop
			}//end of if
629 630 631 632 633 634
			else
			{
				std::cerr << "<<ERROR>> you must provide a vector of parameters or a file to load with the bootstrap" << std::endl;
			}
		}
	}
635
	else
636
	{
637
		std::cout << "<<DEBUG>> No gobal bootstrap option. Bootstraping per function" << std::endl;
638
		for(unsigned int i=0; i<fs.size(); ++i)
639 640 641 642
		{
			fs[i]->bootstrap(d, fs_args[i]);
		}
	}
643
}
644

645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813
void compound_function::setMin(const vec& min) 
{
	function::setMin(min);
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		fs[i]->setMin(min);
	}
}
void compound_function::setMax(const vec& max) 
{
	function::setMax(max);
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		fs[i]->setMax(max);
	}
}

//! Number of parameters to this non-linear function
int compound_function::nbParameters() const
{
	int nb_params = 0;
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		if(!is_fixed[i]) {
			nb_params += fs[i]->nbParameters();
		}
	}
	return nb_params;
}

//! Get the vector of parameters for the function
vec compound_function::parameters() const
{
	vec params(nbParameters());
	int current_i = 0;
	for(unsigned int f=0; f<fs.size(); ++f)
	{
		int f_size = fs[f]->nbParameters();

		// Handle when there is no parameters to include
		if(f_size > 0 && !is_fixed[f])
		{
			vec f_params = fs[f]->parameters();
			for(int i=0; i<f_size; ++i)
			{
				params[i + current_i] = f_params[i];
			}

			current_i += f_size;
		}
	}

	return params;
}

//! Get the vector of min parameters for the function
vec compound_function::getParametersMin() const
{
	vec params(nbParameters());
	int current_i = 0;
	for(unsigned int f=0; f<fs.size(); ++f)
	{
		int f_size = fs[f]->nbParameters();

		// Handle when there is no parameters to include
		if(f_size > 0 && !is_fixed[f])
		{
			vec f_params = fs[f]->getParametersMin();
			for(int i=0; i<f_size; ++i)
			{
				params[i + current_i] = f_params[i];
			}

			current_i += f_size;
		}
	}

	return params;
}

//! Get the vector of min parameters for the function
vec compound_function::getParametersMax() const
{
	vec params(nbParameters());
	int current_i = 0;
	for(unsigned int f=0; f<fs.size(); ++f)
	{
		int f_size = fs[f]->nbParameters();

		// Handle when there is no parameters to include
		if(f_size > 0 && !is_fixed[f])
		{
			vec f_params = fs[f]->getParametersMax();
			for(int i=0; i<f_size; ++i)
			{
				params[i + current_i] = f_params[i];
			}

			current_i += f_size;
		}
	}

	return params;
}

//! Update the vector of parameters for the function
void compound_function::setParameters(const vec& p) 
{
	int current_i = 0;
	for(unsigned int f=0; f<fs.size(); ++f)
	{
		int f_size = fs[f]->nbParameters();

		// Handle when there is no parameters to include
		if(f_size > 0 && !is_fixed[f])
		{
			vec f_params(f_size);
			for(int i=0; i<f_params.size(); ++i)
			{
				f_params[i] = p[i + current_i];
			}

			fs[f]->setParameters(f_params);
			current_i += f_size;
		}
	}
}
		
void compound_function::save_body(std::ostream& out, const arguments& args) const
{
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		fs[i]->save_body(out, args);
		out << std::endl;
	}

	function::save_body(out, args);
}
		
void compound_function::save_call(std::ostream& out, const arguments& args) const
{
	bool is_cpp    = args["export"] == "C++";
	bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
	bool is_matlab = args["export"] == "matlab";

	// This part is export specific. For ALTA, the coefficients are just
	// dumped as is with a #FUNC {plugin_name} header.
	//
	// For C++ export, the function call should be done before hand and
	// the line should look like:
	//   res += call_i(x);
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		if(i != 0 && (is_cpp || is_matlab || is_shader))
		{
			out << "\tres += ";
		}

		fs[i]->save_call(out, args);

		if(is_cpp || is_matlab || is_shader)
		{
			out << ";" << std::endl;
		}
	}

	function::save_call(out, args);
}

814

815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832

// Return the parameters of the product of G1 and G2.
static const parameters
product_parameters(const ptr<nonlinear_function>& g1,
                   const ptr<nonlinear_function>& g2)
{
    parameters result = g1->parametrization();

    if(g1->parametrization().input_parametrization() !=
       g2->parametrization().input_parametrization())
    {
        result = alta::parameters(6, result.dimY(),
                                  params::CARTESIAN,
                                  result.output_parametrization());
    }

    return result;
}
833 834 835

/*--- Product functions implementation ----*/

836
product_function::product_function(const ptr<nonlinear_function>& g1,
837
                                   const ptr<nonlinear_function>& g2,
838 839 840
                                   bool is_g1_fixed, bool is_g2_fixed)
: nonlinear_function(product_parameters(g1, g2)),
  f1( g1 ),
841 842
	f2( g2 ),
	_is_fixed( std::pair<bool,bool>( is_g1_fixed, is_g2_fixed) )
843 844 845
{
}

846 847 848 849 850
product_function::~product_function()
{
}


851 852 853 854 855 856
vec product_function::operator()(const vec& x) const
{
	return value(x);
}
vec product_function::value(const vec& x) const
{
857 858
	// Convert input space to the 1rst function parametrization and compute the
	// output value
859 860 861
	vec xf1(f1->parametrization().dimX());
	params::convert(&x[0], parametrization().input_parametrization(),
                  f1->parametrization().input_parametrization(), &xf1[0]);
862 863 864 865
	vec f1res = f1->value(xf1);

	// Convert input space to the 2nd function parametrization and compute the
	// output value
866 867 868
	vec xf2(f2->parametrization().dimX());
	params::convert(&x[0], parametrization().input_parametrization(),
                  f2->parametrization().input_parametrization(), &xf2[0]);
869 870
	vec f2res = f2->value(xf2);

871 872 873 874 875
	vec res = product(f1res, f2res);

	//std::cout << f1res << " :::::::::: " << res << std::endl;
	//std::cout << f2res << std::endl;
	return res;
876 877 878 879
}
		
bool product_function::load(std::istream& in)
{
880
	bool loaded_f1 = false,loaded_f2 = false;
881 882 883
	std::streampos pos = in.tellg();

	// Load the first function
884
	if(f1) {
885 886 887 888 889 890 891 892 893
		loaded_f1 = f1->load(in);
		if(! loaded_f1)
		{
			in.seekg(pos);
		}
	}

	pos = in.tellg();
	// Load the second function
894
	if(f2) {
895 896 897 898 899 900 901 902 903 904
		loaded_f2 = f2->load(in);
		if(! loaded_f2)
		{
			in.seekg(pos);
		}
	}

	return loaded_f1 || loaded_f2;
}

905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949
void product_function::save_body(std::ostream& out, const arguments& args) const
{
	f1->save_body(out, args);
	out << std::endl;
	
	f2->save_body(out, args);
	out << std::endl;

	function::save_body(out, args);
}

void product_function::save_call(std::ostream& out, const arguments& args) const
{
	bool is_cpp    = args["export"] == "C++";
	bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
	bool is_matlab = args["export"] == "matlab";

	// This part is export specific. For ALTA, the coefficients are just
	// dumped as is with a #FUNC {plugin_name} header.
	//
	// For C++ export, the function call should be done before hand and
	// the line should look like:
	//   res += call_i(x);
	if(is_cpp || is_matlab || is_shader)
	{
		out << "(";
	}

	f1->save_call(out, args);

	if(is_cpp || is_matlab || is_shader)
	{
		out << " * ";
	}

	f2->save_call(out, args);

	if(is_cpp || is_matlab || is_shader)
	{
		out << ")" << std::endl;
	}

	function::save_call(out, args);
}

950
void product_function::bootstrap(const ptr<data> d, const arguments& args)
951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985
{
	const bool global_bootstrap = args.is_defined("bootstrap");

	// Check if the bootstrap is global
	if(global_bootstrap)
	{
		if(args.is_vec("bootstrap"))
		{
			vec p = args.get_vec("bootstrap", nbParameters());
			setParameters(p);
		}
		else
		{
			std::ifstream file;
			file.open(args["bootstrap"].c_str()) ;
			if(file.is_open())
			{
				// Skip the header
				std::string line ;
				do
				{
					std::getline(file, line) ;
				}
				while(line != "#ALTA HEADER END");

				std::streampos pos = file.tellg();

				// If the first function cannot be loaded, put the input stream
				// in the previous state and bootstrap normaly this function.
				if(!f1->load(file))
				{
					file.seekg(pos);

					// Bootstrap the function as if it was not loaded
					f1->bootstrap(d, args);
986
					std::cout << "<<DEBUG>> Unable to load first function of product, regular bootstraping" << std::endl;
987 988 989 990 991 992 993 994 995 996
				}
				
				// If the second function cannot be loaded, put the input stream
				// in the previous state and bootstrap normaly this function.
				if(!f2->load(file))
				{
					file.seekg(pos);

					// Bootstrap the function as if it was not loaded
					f2->bootstrap(d, args);
997
					std::cout << "<<DEBUG>> Unable to load second function of product, regular bootstraping" << std::endl;
998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014
				}
			}
			else
			{
				std::cerr << "<<ERROR>> you must provide a vector of parameters or a file to load with the bootstrap" << std::endl;
			}
		}
	}
	else
	{
		f1->bootstrap(d, args);
		f2->bootstrap(d, args);
	}
}

void product_function::setMin(const vec& min) 
{
1015
	function::setMin(min);
1016 1017 1018 1019 1020 1021
	f1->setMin(min);
	f2->setMin(min);
}

void product_function::setMax(const vec& max) 
{
1022
	function::setMax(max);
1023 1024 1025 1026 1027 1028
	f1->setMax(max);
	f2->setMax(max);
}

int product_function::nbParameters() const
{
1029 1030 1031 1032 1033 1034 1035
	int num_parameters = 0;
	
	num_parameters +=	(! _is_fixed.first) ? (f1->nbParameters()) : (0) ;

	num_parameters +=	(! _is_fixed.second) ? (f2->nbParameters()) : (0) ;

	return num_parameters;
1036 1037 1038 1039
}

vec product_function::parameters() const
{
1040 1041 1042 1043
	// int nb_f1_params = f1->nbParameters();
	// int nb_f2_params = f2->nbParameters();
	// int nb_params = nb_f1_params + nb_f2_params;
	int const nb_params = product_function::nbParameters();	
1044

1045 1046
	int nb_f1_params = 0;
	
1047
	vec params(nb_params);
1048
	if( !_is_fixed.first )
1049
	{
1050 1051 1052 1053 1054 1055 1056 1057 1058
		//F1 is active save the number of paramters of f1 
		nb_f1_params = f1->nbParameters();

		vec f1_params = f1->parameters();
		for(int i=0; i<nb_f1_params; ++i)
		{
			params[i] = f1_params[i];
		}

1059 1060
	}

1061
	if( ! _is_fixed.second )
1062
	{
1063 1064 1065 1066 1067
		int const nb_f2_params = f2->nbParameters();
		vec f2_params = f2->parameters();
		for(int i=0; i<nb_f2_params; ++i)
		{
			params[i+nb_f1_params] = f2_params[i];
1068
		}
1069 1070 1071 1072
	}

	return params;
}
1073 1074 1075 1076

void product_function::setParameters(const vec& p)
{

1077 1078 1079
	int nb_f1_params = 0;

	if( ! _is_fixed.first )
1080
	{
1081 1082 1083 1084 1085 1086 1087 1088
		nb_f1_params = f1->nbParameters();
	
		vec f1_params(nb_f1_params);
		for(int i=0; i<nb_f1_params; ++i)
		{
			f1_params[i] = p[i];
		}

1089
		f1->setParameters(f1_params);	    
1090 1091
	}

1092
	if( ! _is_fixed.second )
1093
	{
1094 1095 1096 1097 1098 1099 1100 1101
		int const nb_f2_params = f2->nbParameters();
		vec f2_params(nb_f2_params);
		for(int i=0; i<nb_f2_params; ++i)
		{
			f2_params[i] = p[i+nb_f1_params];
		}
		
		f2->setParameters(f2_params);			
1102
	}
1103 1104


1105 1106
}

1107
//! Get the vector of max parameters for the function
1108 1109
vec product_function::getParametersMax() const
{
1110 1111 1112
	int const nb_params = product_function::nbParameters();	

	int nb_f1_params = 0;
1113 1114 1115

	vec params(nb_params);

1116
	if( !_is_fixed.first )
1117
	{
1118 1119 1120 1121 1122 1123
		nb_f1_params = f1->nbParameters();
		vec f1_params = f1->getParametersMax();
		for(int i=0; i<nb_f1_params; ++i)
		{
			params[i] = f1_params[i];
		}
1124 1125
	}

1126
	if( ! _is_fixed.second )
1127
	{
1128 1129 1130 1131 1132 1133 1134
		int const nb_f2_params = f2->nbParameters();
		vec f2_params = f2->getParametersMax();
		for(int i=0; i<nb_f2_params; ++i)
		{
			params[i+nb_f1_params] = f2_params[i];
		}

1135 1136 1137 1138 1139 1140 1141 1142
	}

	return params;
}

//! Get the vector of min parameters for the f1tion
vec product_function::getParametersMin() const
{
1143 1144 1145
	int nb_f1_params = 0;

	int const nb_params = product_function::nbParameters();	
1146 1147 1148

	vec params(nb_params);

1149
	if( ! _is_fixed.first )
1150
	{
1151 1152 1153 1154 1155 1156 1157
		nb_f1_params = f1->nbParameters();

		vec f1_params = f1->getParametersMin();
		for(int i=0; i<nb_f1_params; ++i)
		{
			params[i] = f1_params[i];
		}
1158 1159
	}

1160
	if( ! _is_fixed.second )
1161
	{
1162 1163 1164 1165 1166 1167
		int const nb_f2_params = f2->nbParameters();
		vec f2_params = f2->getParametersMin();
		for(int i=0; i<nb_f2_params; ++i)
		{
			params[i+nb_f1_params] = f2_params[i];
		}	
1168 1169
	}

1170 1171 1172

	

1173 1174 1175
	return params;
}

1176 1177
vec product_function::parametersJacobian(const vec& x) const
{
1178 1179 1180 1181 1182 1183 1184 1185
//ß	std::cout << "ENTERING parametersJacobian " __FILE__ << " " << __LINE__ << std::endl;


	int const nb_f1_params = f1->nbParameters();
	int const nb_f2_params = f2->nbParameters();
	//int nb_params = nb_f1_params + nb_f2_params;
	int const nb_params = product_function::nbParameters();

1186 1187

	// Convert the input value x to the input space of the f1tion
1188 1189 1190 1191 1192 1193
	vec xf2(f2->parametrization().dimX());
	params::convert(&x[0], parametrization().input_parametrization(),
                  f2->parametrization().input_parametrization(), &xf2[0]);
	vec xf1(f1->parametrization().dimX());
	params::convert(&x[0], parametrization().input_parametrization(),
                  f1->parametrization().input_parametrization(), &xf1[0]);
1194

1195
	//Value of each function for the given x
1196 1197
	vec f1_value = f1->value(xf1);
	vec f2_value = f2->value(xf2);
1198

1199 1200 1201 1202
	//Jacobien of each function for the given x
	vec f1_jacobian = f1->parametersJacobian(xf1);
	vec f2_jacobian = f2->parametersJacobian(xf2);

1203 1204
	// F = f2nel; f = f1tion
	// d(F * f)(x) /dp = F(x) df(x) /dp + f(x) dF(x) / dp
1205
	//vec jac(nb_params*parametrization().dimY());
1206 1207
	
	//Forcing Zero by default
1208
	vec jac= vec::Zero( nb_params * parametrization().dimY() );
1209 1210


1211
	for(int y=0; y<parametrization().dimY(); ++y)
1212
	{
1213
		if( ! _is_fixed.first)
1214
		{
1215 1216 1217 1218
			for(int i=0; i<nb_f1_params; ++i)
			{
				jac[y*nb_params + i] = f1_jacobian[y*nb_f1_params + i] * f2_value[y];
			}	
1219
		}
1220 1221
		
		if( ! _is_fixed.second )
1222
		{
1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234
			for(int i=0; i<nb_f2_params; ++i)
			{
				if( _is_fixed.first )
				{
					jac[y*nb_params + i ] = f2_jacobian[y*nb_f2_params + i] * f1_value[y];
				}
				else
				{
					jac[y*nb_params + (i+nb_f1_params)] = f2_jacobian[y*nb_f2_params + i] * f1_value[y];	
				}
				
			}				
1235 1236 1237 1238 1239 1240
		}
	}

	return jac;
}

1241 1242
nonlinear_function* product_function::first() const
{
1243
	return f1.get();
1244 1245 1246 1247
}

nonlinear_function* product_function::second() const
{
1248
	return f2.get();
1249
}