function.cpp 27.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

   Copyright (C) 2015 CNRS
   Copyright (C) 2013, 2014 Inria

   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
/*--- Functions implementation ----*/

20
void function::bootstrap(const ptr<data>, const arguments& args)
21 22 23 24 25 26 27
{
    // 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())
        {
28 29 30 31
            if(!load(in))
				{
					std::cerr << "<<ERROR>> cannot bootstrap from file \"" << args["bootstrap"] << "\"" << std::endl;
				}
32 33 34 35
        }
    }
}

36 37
void function::save(const std::string& filename, const arguments& args) const
{
PACANOWSKI Romain's avatar
PACANOWSKI Romain committed
38 39 40 41 42
	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";
43

44
	// Open the file
Laurent Belcour's avatar
Laurent Belcour committed
45
	std::ofstream file(filename.c_str());
46 47 48 49 50
	if(!file.is_open())
	{
		std::cerr << "<<ERROR>> unable to open output file for writing" << std::endl;
	}

51 52 53 54 55 56 57 58 59 60 61 62
	// 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;
	}
63

64
	// Save common header
65
	save_header(file, args);
66 67

	// Save function definition
68
	save_body(file, args);
69

70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
	if(is_cpp)
	{
		file << "vec brdf(const vec& in, const vec& file)" << std::endl;
		file << "{" << std::endl;
		file << "\tvec res(" << dimY() << ");" << std::endl;
		file << "\t";
	}
	else if(is_matlab)
	{
		file << "function res = brdf(in, file)" << std::endl;
		file << "{" << std::endl;
		file << "\tres = zeros(" << dimY() << ");" << std::endl;
		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 = ";
	}
90 91

	// Save fit data
92 93 94 95 96 97 98 99 100 101 102 103 104 105
	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;
106

107 108 109 110
	if(is_explorer)
	{
		file << "::end shader" << std::endl;
	}
111 112 113 114 115 116 117 118 119
}
		
//! \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"))
	{
120
		out << "#ALTA FUNC HEADER" << std::endl;
121 122 123
		out << "#CMD " << args.get_cmd() << std::endl;
		out << "#DIM " << _nX << " " << _nY << std::endl;
		out << "#PARAM_IN  " << params::get_name(input_parametrization()) << std::endl;
124
		//out << "#PARAM_OUT " << params::get_name(output_parametrization()) << std::endl;*
125 126 127 128
		if(args.is_defined("export-append")) 
		{
			out << args["export-append"] << std::endl;
		}
129
		out << "#ALTA HEADER END" << std::endl;
130 131 132 133 134 135 136
		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.
137
void function::save_body(std::ostream&, const arguments&) const
138 139
{

140

141 142 143 144 145
}

//! \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.
146
void function::save_call(std::ostream&, const arguments&) const
147 148 149 150
{

}

151
//! \brief L2 norm to data.
152
double function::L2_distance(const ptr<data>& d) const
153
{
154
	double l2_dist = 0.0;
155 156 157
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
158
		vec x(dimX()), y(dimY());
159

160 161 162 163 164 165 166 167 168 169
		if(input_parametrization() == params::UNKNOWN_INPUT)
		{
			memcpy(&x[0], &dat[0], dimX()*sizeof(double));
		}
		else
		{
			params::convert(&dat[0], d->input_parametrization(), input_parametrization(), &x[0]);
		}
		memcpy(&y[0], &dat[d->dimX()], dimY()*sizeof(double));

170
		l2_dist += std::pow(norm(y-value(x)), 2);
171 172
	}
	l2_dist = std::sqrt(l2_dist / d->size());
173
	return l2_dist;
174 175 176
}

//! \brief Linf norm to data.
177
double function::Linf_distance(const ptr<data>& d) const
178
{
179 180
	vec mean = vec::Zero(dimY());
	vec var  = vec::Zero(dimY());
181

182
#ifdef DEBUG
183
	std::cout << "<<DEBUG>> input param here = " << params::get_name(input_parametrization()) << std::endl;
184
#endif
185

186
	double linf_dist = 0.0;
187 188 189
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
190
		vec x(dimX()), y(dimY());
191

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
        // Convert the position of the data sample to the parametrization
        // of the function.
        if(input_parametrization() == params::UNKNOWN_INPUT)
        {
            memcpy(&x[0], &dat[0], dimX()*sizeof(double));
        }
        else
        {
            params::convert(&dat[0], d->input_parametrization(), input_parametrization(), &x[0]);
        }

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

		// Take the componentwise-max of the two vectors.
		const vec v = value(x);
        for(int j=0; j<d->dimY(); ++j)
210
		{
211
			linf_dist = std::max<double>(linf_dist, std::abs(y[j]-v[j]));
212
		}
213

214 215
		// Compute the mean
		mean += (y-v) / double(d->size());
216 217
	}
	
218
	// Compute the standard deviation with respect to the mean error
219 220 221 222
	for(int i=0; i<d->size(); ++i)
	{
		vec dat = d->get(i);
		vec x(dimX()), y(d->dimY()), val(dimY());
223 224 225 226 227 228 229 230 231 232 233
		
        // Convert the position of the data sample to the parametrization
        // of the function.
        if(input_parametrization() == params::UNKNOWN_INPUT)
        {
            memcpy(&x[0], &dat[0], dimX()*sizeof(double));
        }
        else
        {
            params::convert(&dat[0], d->input_parametrization(), input_parametrization(), &x[0]);
        }
234

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

239
		val = value(x);
240 241 242 243 244
		for(int j=0; j<d->dimY(); ++j) 
		{ 
			y[j] = dat[d->dimX()+j]; 
			var[j] += pow(mean[j] - (val[j]-y[j]), 2) / double(d->size());
		}
245 246
	}

247 248
	std::cout << "<<INFO>> Mean = " << mean << ", Var = " << var << std::endl;

249
	return linf_dist;
250
}
251 252


253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 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

/*--- Non-linear functions implementation ----*/
		
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);
}

320
void nonlinear_function::bootstrap(const ptr<data> d, const arguments& args)
321 322 323 324 325 326 327 328 329 330 331
{
    if(args.is_vec("bootstrap"))
    {
        vec p = args.get_vec("bootstrap", nbParameters());
        setParameters(p);
    }
    else
    {
        function::bootstrap(d, args);
    }
}
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351
		
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;
}
352

353 354 355 356 357 358 359 360 361 362


/*--- Compound functions implementation ----*/

vec compound_function::operator()(const vec& x) const
{
	return value(x);
}
vec compound_function::value(const vec& x) const
{
363
	vec res = vec::Zero(dimY());
364 365
	for(unsigned int i=0; i<fs.size(); ++i)
	{
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387
		vec temp_x(fs[i]->dimX());
		params::convert(&x[0], input_parametrization(), fs[i]->input_parametrization(), &temp_x[0]);
		res = res + fs[i]->value(temp_x);
	}
	return res;
}

vec compound_function::parametersJacobian(const vec& x) const
{
	int nb_params = nbParameters();
	vec jac = vec::Zero(nb_params*dimY());

	int start_i = 0;

	// Export the sub-Jacobian for each function
	for(unsigned int f=0; f<fs.size(); ++f)
	{
		nonlinear_function* func = fs[f];
		int nb_f_params = func->nbParameters(); 

		// Only export Jacobian if there are non-linear parameters
		if(nb_f_params > 0 && !is_fixed[f])
388
		{
389 390 391 392 393 394 395 396 397 398 399 400 401 402

			vec temp_x(func->dimX());
			params::convert(&x[0], input_parametrization(), func->input_parametrization(), &temp_x[0]);
			vec func_jac = func->parametersJacobian(temp_x);

			for(int i=0; i<nb_f_params; ++i)
			{
				for(int y=0; y<dimY(); ++y)
				{
					jac[y*nb_params + (i+start_i)] = func_jac[y*nb_f_params + i];
				}
			}

			start_i += nb_f_params;
403 404
		}
	}
405 406

	return jac;
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
}

void compound_function::push_back(nonlinear_function* f, const arguments& f_args)
{
	// Update the input param
	if(input_parametrization() == params::UNKNOWN_INPUT)
	{
		setParametrization(f->input_parametrization());
	}
	else if(input_parametrization() != f->input_parametrization())
	{
		setParametrization(params::CARTESIAN);
	}

	// Update the output param
	if(output_parametrization() == params::UNKNOWN_OUTPUT)
	{
		setParametrization(f->output_parametrization());
	}
	else if(output_parametrization() != f->output_parametrization())
	{
		std::cerr << "Creating a compound function with different output dimensions, this is not allowed" << std::endl;
		throw;
	}

	fs.push_back(f);
	fs_args.push_back(f_args);
	is_fixed.push_back(f_args.is_defined("fixed"));
}

nonlinear_function* compound_function::operator[](int i) const
{
#ifdef DEBUG
	assert(i >= 0 && i < fs.size());
#endif
	return fs[i];
}
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
		
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++;
		}
	}
465
	std::cout << "<<DEBUG>> number of correct function loaded in compound: " << nb_good  << " / " << fs.size() << std::endl;
466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500
	return nb_good > 0;
}

void compound_function::setParametrization(params::input new_param)
{
	if(new_param == params::UNKNOWN_INPUT)
		return;

	// If there is more than one parametrization defined to this function, convert
	// it to the parametrization that conserves the most degrees of freedom. Right now
	// this is the CARTESIAN param. It might change in future release.
	if(input_parametrization() != new_param)
	{
		function::setParametrization(params::CARTESIAN);
		function::setDimX(6);
	}

	for(unsigned int i=0; i<fs.size(); ++i)
	{
		if(fs[i]->input_parametrization() == params::UNKNOWN_INPUT)
		{
			fs[i]->setParametrization(new_param);
			fs[i]->setDimX(params::dimension(new_param));
		}
	}
}
		
void compound_function::setParametrization(params::output new_param)
{
	parametrized::setParametrization(new_param);
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		fs[i]->setParametrization(new_param);
	}
}
501

502
void compound_function::bootstrap(const ::ptr<data> d, const arguments& args)
503
{
504
	const bool global_bootstrap = args.is_defined("bootstrap");
505

506 507 508 509 510 511
	// Check if the bootstrap is global
	if(global_bootstrap)
	{
		if(args.is_vec("bootstrap"))
		{
			vec p = args.get_vec("bootstrap", nbParameters());
512
			std::cout << "<<INFO>> Will use " << p << " as a bootstrap for the non-linear function" << std::endl;
513
			setParameters(p);
514 515 516 517 518 519 520 521
		
			for(unsigned int i=0; i<fs.size(); ++i)
			{
				if(fs[i]->nbParameters() == 0)
				{
					fs[i]->bootstrap(d, fs_args[i]);
				}
			}
522 523 524 525 526 527 528 529 530 531 532 533 534 535
		}
		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");
536

537 538 539 540 541 542 543
				// 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();
544
					
545 546
					if(dynamic_cast<product_function*>(fs[i]) != NULL)
					{
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571
						product_function* p = dynamic_cast<product_function*>(fs[i]);
						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;
						}
572
					}
573 574 575 576 577 578 579
					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);
580

581 582 583 584
							// 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;
						}
585 586 587 588 589 590 591 592 593
					}
				}
			}
			else
			{
				std::cerr << "<<ERROR>> you must provide a vector of parameters or a file to load with the bootstrap" << std::endl;
			}
		}
	}
594
	else
595
	{
596
		for(unsigned int i=0; i<fs.size(); ++i)
597 598 599 600
		{
			fs[i]->bootstrap(d, fs_args[i]);
		}
	}
601
}
602

603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 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
void compound_function::setDimX(int nX) 
{
	if(input_parametrization() == params::UNKNOWN_INPUT)
	{
		function::setDimX(nX);
	}
	
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		if(fs[i]->input_parametrization() == params::UNKNOWN_INPUT)
		{
			fs[i]->setDimX(nX);
		}
	}
}
		
void compound_function::setDimY(int nY)
{
	function::setDimY(nY);
	for(unsigned int i=0; i<fs.size(); ++i)
	{
		fs[i]->setDimY(nY);
	}
}

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);
}

797 798 799 800 801 802 803 804


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

product_function::product_function(nonlinear_function* g1, nonlinear_function* g2) 
{
	this->f1 = g1;
	this->f2 = g2;
805 806 807 808 809 810 811 812 813 814

	// If the two parametrization are different, use the CARTESIAN parametrization
	// as the input parametrization, then do the convertion for all the functions.
	if(g1->input_parametrization() != g2->input_parametrization())
	{
		function::setParametrization(params::CARTESIAN);
		function::setDimX(6);
	}
	else
	{
815
		function::setParametrization(g1->input_parametrization());
816 817
		function::setDimX(g1->dimX());
	}
818 819 820 821 822 823 824 825
}

vec product_function::operator()(const vec& x) const
{
	return value(x);
}
vec product_function::value(const vec& x) const
{
826 827 828 829 830 831 832 833
	// Convert input space to the 1rst function parametrization and compute the
	// output value
	vec xf1(f1->dimX());
	params::convert(&x[0], input_parametrization(), f1->input_parametrization(), &xf1[0]);
	vec f1res = f1->value(xf1);

	// Convert input space to the 2nd function parametrization and compute the
	// output value
834
	vec xf2(f2->dimX());
835
	params::convert(&x[0], input_parametrization(), f2->input_parametrization(), &xf2[0]);
836 837
	vec f2res = f2->value(xf2);

838 839 840 841 842
	vec res = product(f1res, f2res);

	//std::cout << f1res << " :::::::::: " << res << std::endl;
	//std::cout << f2res << std::endl;
	return res;
843 844 845 846
}
		
bool product_function::load(std::istream& in)
{
847
	bool loaded_f1 = false,loaded_f2 = false;
848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873
	std::streampos pos = in.tellg();

	// Load the first function
	if(f1 != NULL)
	{
		loaded_f1 = f1->load(in);
		if(! loaded_f1)
		{
			in.seekg(pos);
		}
	}

	pos = in.tellg();
	// Load the second function
	if(f2 != NULL)
	{
		loaded_f2 = f2->load(in);
		if(! loaded_f2)
		{
			in.seekg(pos);
		}
	}

	return loaded_f1 || loaded_f2;
}

874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918
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);
}

919
void product_function::bootstrap(const ptr<data> d, const arguments& args)
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 950 951 952 953 954
{
	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);
955
					std::cout << "<<DEBUG>> Unable to load first function of product, regular bootstraping" << std::endl;
956 957 958 959 960 961 962 963 964 965
				}
				
				// 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);
966
					std::cout << "<<DEBUG>> Unable to load second function of product, regular bootstraping" << std::endl;
967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985
				}
			}
			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::setDimX(int nX) 
{
	f1->setDimX(nX);
	f2->setDimX(nX);
986 987
	
	function::setDimX(f1->dimX());
988 989 990 991 992 993
}

void product_function::setDimY(int nY)
{
	f1->setDimY(nY);
	f2->setDimY(nY);
994 995
	
	function::setDimY(f1->dimY());
996 997 998 999
}

void product_function::setMin(const vec& min) 
{
1000
	function::setMin(min);
1001 1002 1003 1004 1005 1006
	f1->setMin(min);
	f2->setMin(min);
}

void product_function::setMax(const vec& max) 
{
1007
	function::setMax(max);
1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030
	f1->setMax(max);
	f2->setMax(max);
}

int product_function::nbParameters() const
{
	return f1->nbParameters() + f2->nbParameters();
}

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

	vec params(nb_params);

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

1031
	vec f2_params = f2->parameters();
1032 1033 1034 1035 1036 1037 1038
	for(int i=0; i<nb_f2_params; ++i)
	{
		params[i+nb_f1_params] = f2_params[i];
	}

	return params;
}
Laurent Belcour's avatar
Laurent Belcour committed
1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054

void product_function::setParameters(const vec& p)
{
	int nb_f1_params = f1->nbParameters();
	int nb_f2_params = f2->nbParameters();

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

	vec f2_params(nb_f2_params);
	for(int i=0; i<nb_f2_params; ++i)
	{
1055
		f2_params[i] = p[i+nb_f1_params];
Laurent Belcour's avatar
Laurent Belcour committed
1056 1057 1058 1059
	}
	f2->setParameters(f2_params);
}

1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107
//! Get the vector of min parameters for the function
vec product_function::getParametersMax() const
{
	int nb_f1_params = f1->nbParameters();
	int nb_f2_params = f2->nbParameters();
	int nb_params = nb_f1_params + nb_f2_params;

	vec params(nb_params);

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

	vec f2_params = f2->getParametersMax();
	for(int i=0; i<nb_f2_params; ++i)
	{
		params[i+nb_f1_params] = f2_params[i];
	}

	return params;
}

//! Get the vector of min parameters for the f1tion
vec product_function::getParametersMin() const
{
	int nb_f1_params = f1->nbParameters();
	int nb_f2_params = f2->nbParameters();
	int nb_params = nb_f1_params + nb_f2_params;

	vec params(nb_params);

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

	vec f2_params = f2->getParametersMin();
	for(int i=0; i<nb_f2_params; ++i)
	{
		params[i+nb_f1_params] = f2_params[i];
	}

	return params;
}

Laurent Belcour's avatar
Laurent Belcour committed
1108 1109 1110 1111 1112 1113 1114
vec product_function::parametersJacobian(const vec& x) const
{
	int nb_f1_params = f1->nbParameters();
	int nb_f2_params = f2->nbParameters();
	int nb_params = nb_f1_params + nb_f2_params;

	// Convert the input value x to the input space of the f1tion
1115 1116 1117 1118
	vec xf2(f2->dimX());
	params::convert(&x[0], input_parametrization(), f2->input_parametrization(), &xf2[0]);
	vec xf1(f1->dimX());
	params::convert(&x[0], input_parametrization(), f1->input_parametrization(), &xf1[0]);
Laurent Belcour's avatar
Laurent Belcour committed
1119

1120 1121
	vec f1_jacobian = f1->parametersJacobian(xf1);
	vec f2_jacobian = f2->parametersJacobian(xf2);
Laurent Belcour's avatar
Laurent Belcour committed
1122

1123 1124
	vec f1_value = f1->value(xf1);
	vec f2_value = f2->value(xf2);
Laurent Belcour's avatar
Laurent Belcour committed
1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149

	// F = f2nel; f = f1tion
	// d(F * f)(x) /dp = F(x) df(x) /dp + f(x) dF(x) / dp
	vec jac(nb_params*_nY);
	for(int y=0; y<_nY; ++y)
	{
		for(int i=0; i<nb_f1_params; ++i)
		{
			jac[y*nb_params + i] = f1_jacobian[y*nb_f1_params + i] * f2_value[y];
		}

		for(int i=0; i<nb_f2_params; ++i)
		{
			jac[y*nb_params + (i+nb_f1_params)] = f2_jacobian[y*nb_f2_params + i] * f1_value[y];
		}
	}

	return jac;
}

//! \brief provide the outout parametrization of the object.
params::output product_function::output_parametrization() const
{
	return f1->output_parametrization();
}
1150 1151 1152
		
void product_function::setParametrization(params::input new_param)
{
1153 1154 1155 1156 1157 1158 1159 1160 1161 1162
	if(f1->input_parametrization() == params::UNKNOWN_INPUT) 
	{
		f1->setParametrization(new_param); 
		f1->setDimX(params::dimension(new_param));
	}
	if(f2->input_parametrization() == params::UNKNOWN_INPUT) 
	{
		f2->setParametrization(new_param);
		f2->setDimX(params::dimension(new_param));
	}
1163 1164 1165 1166 1167 1168 1169
}

void product_function::setParametrization(params::output new_param)
{
	f1->setParametrization(new_param);
	f2->setParametrization(new_param);
}
1170 1171 1172 1173 1174 1175 1176 1177 1178 1179

nonlinear_function* product_function::first() const
{
	return f1;
}

nonlinear_function* product_function::second() const
{
	return f2;
}