function.cpp 5.83 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
#include "function.h"

#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>

#include <core/common.h>

ALTA_DLL_EXPORT function* provide_function()
{
    return new abc_function();
}

// Overload the function operator
vec abc_function::operator()(const vec& x) const 
{
	return value(x);
}
vec abc_function::value(const vec& x) const 
{

	// Compute the Shadow term to init res
	vec res = vec::Zero(dimY());
	const double hn = 1.0 - x[0];

	for(int i=0; i<dimY(); ++i)
	{
32 33 34 35 36 37 38 39
		if(hn > 0.0 && hn < 1.0)
		{
			res[i] = _a[i] / pow(1.0 + _b[i]*hn, _c[i]);
		}
		else
		{
			res[i] = 0.0;
		}
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
	}
	return res;
}

// Reset the output dimension
void abc_function::setDimY(int nY)
{
    _nY = nY ;

    // Update the length of the vectors
    _a.resize(_nY) ;
    _b.resize(_nY) ;
    _c.resize(_nY) ;
}

//! Number of parameters to this non-linear function
int abc_function::nbParameters() const 
{
	return 3*dimY();
}

//! Get the vector of parameters for the function
vec abc_function::parameters() const 
{
	vec res(nbParameters());
	for(int i=0; i<dimY(); ++i)
	{
		res[i*3 + 0] = _a[i];
		res[i*3 + 1] = _b[i];
		res[i*3 + 2] = _c[i];
	}
	return res;
}

//! \brief get the min values for the parameters
vec abc_function::getParametersMin() const
{
	vec res(nbParameters());
	for(int i=0; i<dimY(); ++i)
	{
		res[i*3 + 0] = 0.0;
		res[i*3 + 1] = 0.0;
		res[i*3 + 2] = 0.0;
	}

	return res;
}

//! Update the vector of parameters for the function
void abc_function::setParameters(const vec& p) 
{
	for(int i=0; i<dimY(); ++i)
	{
		_a[i]  = p[i*3 + 0];
		_b[i]  = p[i*3 + 1];
		_c[i]  = p[i*3 + 2];
	}
}

//! Obtain the derivatives of the function with respect to the 
//! parameters
//! \todo finish. 
vec abc_function::parametersJacobian(const vec& x) const 
{
    vec jac(dimY()*nbParameters());
	 for(int i=0; i<dimY(); ++i)
	 {
		 for(int j=0; j<dimY(); ++j)
		 {
			 if(i == j)
			 {
				 const double hn = 1.0 - x[0];
				 const double f  = 1.0 + _b[i]*hn;
				 const double denom = pow(f, _c[i]);
				 const double fact  = 1.0 / denom;
				 const double fact2 = fact*fact;

				 // df / da
				 jac[i*nbParameters() + j*3+0] = fact;

				 // df / db
				 jac[i*nbParameters() + j*3+1] = - _a[i] * hn *_c[i] * pow(f, _c[i]-1.0) * fact2;

				 // df / dc
				 if(f > 0.0)
				 {
					 jac[i*nbParameters() + j*3+2] = - _a[i] * log(f) * fact;
				 }
				 else
				 {
					jac[i*nbParameters() + j*3+2] = 0.0;
				 }
			 }
			 else
			 {
				 jac[i*nbParameters() + j*3+0] = 0.0;
				 jac[i*nbParameters() + j*3+1] = 0.0;
				 jac[i*nbParameters() + j*3+2] = 0.0;
			 }
		 }
	 }

    return jac;
}
		
void abc_function::bootstrap(const data* d, const arguments& args)
{
	for(int i=0; i<dimY(); ++i)
	{
		_a[i]  = 1.0;
		_b[i]  = 1.0;
		_c[i]  = 1.0;
	}
153 154 155 156 157 158 159 160 161 162 163 164 165
	
	if(args.is_defined("param"))
	{
		setParametrization(params::parse_input(args["param"]));
	}
	else
	{
		setParametrization(params::COS_TH);
	}
	if(params::dimension(input_parametrization()) != 1)
	{
		std::cerr << "<<ERROR>> the parametrization specifed in the file for the ABC model is incorrect" << std::endl;
	}
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
}

//! Load function specific files
bool abc_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_lafortune
	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; 
        return false;
	}

	in >> token;
   if(token.compare("nonlinear_function_abc") != 0) 
	{
		std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl; 
        return false;
	}
198 199 200 201 202 203 204 205 206 207 208 209 210
	
	in >> token;
	if(token.compare("#PARAM") != 0) 
	{ 
		std::cerr << "<<ERROR>> parsing the stream. The #PARAM is not the next line defined." << std::endl; 
        return false;
	}
	in >> token;
	setParametrization(params::parse_input(token));
	if(params::dimension(input_parametrization()) != 1)
	{
		std::cerr << "<<ERROR>> the parametrization specifed in the file for the ABC model is incorrect" << std::endl;
	}
211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230

	// Parse the lobe
	for(int i=0; i<_nY; ++i)
	{

		in >> token >> _a[i];
		in >> token >> _b[i];
		in >> token >> _c[i];
	}
    return true;
}


void abc_function::save_call(std::ostream& out, const arguments& args) const
{
    bool is_alta   = !args.is_defined("export") || args["export"] == "alta";

    if(is_alta)
    {
		out << "#FUNC nonlinear_function_abc" << std::endl ;
231
		out << "#PARAM " << params::get_name(input_parametrization()) << std::endl;
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 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

		 for(int i=0; i<_nY; ++i)
		 {
			 out << "a  " << _a[i]  << std::endl;
			 out << "b  " << _b[i]  << std::endl;
			 out << "c  " << _c[i]  << std::endl;
		 }
	
		 out << std::endl;
	 }
	 else
	 {
		 out << "abc(L, V, N, X, Y, vec3(";
		 for(int i=0; i<_nY; ++i)
		 {
			 out << _a[i];
			 if(i < _nY-1) { out << ", "; }
		 }

		 out << "), vec3(";
		 for(int i=0; i<_nY; ++i)
		 {
			 out << _b[i];
			 if(i < _nY-1) { out << ", "; }
		 }

		 out << "), vec3(";
		 for(int i=0; i<_nY; ++i)
		 {
			 out << _c[i];
			 if(i < _nY-1) { out << ", "; }
		 }
		 out << "))";
	 }

}

void abc_function::save_body(std::ostream& out, const arguments& args) const
{
    bool is_shader = args["export"] == "shader" || args["export"] == "explorer";

    if(is_shader)
    {
        out << "vec3 abc(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y, vec3 a, vec3 b, vec3 c)" << std::endl;
        out << "{" << std::endl;
        out << "\tvec3  H   = normalize(L + V);" << std::endl;
        out << "\tfloat hn  = dot(H,N);" << std::endl;
		  out << "\t" << std::endl;
        out << "\treturn a / pow(vec3(1.0f) + b*hn, c);" << std::endl;
        out << "}" << std::endl;
    }
}