text:
#include <iostream>
#include <fstream>
#include <ctime>
#include <boost/random/mersenne_twister.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>

#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/uniform_real_distribution.hpp>
#include <sys/time.h>

#include <vector>
#include <string>
#include <algorithm>
#include <numeric>
#include <queue>

#include <math.h>
#include <matrix.h>
#include <mex.h>

boost::random::mt19937 gen;

using namespace std;

int collisions (int &dimx_xp, deque<int> &xist_gene, deque<int> &tsix_gene,
        int &sum_xp, int &sum_tp, boost::random::uniform_int_distribution<> &dist)
{
    int q;
    for (q=0; q<dimx_xp; q++)
            {
            if ((xist_gene[q] + tsix_gene[q])>1) {
                if (dist(gen)==0) {
                    xist_gene[q]=0;
                    sum_xp--;
                } else {
                    tsix_gene[q]=0;
                    sum_tp--;
                }
            }
        }
}

int elongation (int &dimx_xp, deque<int> &xist_gene, deque<int> &tsix_gene, int &sum_xp, int &sum_tp,
        bool &xpol_ini, bool &tpol_ini, int &xist_rna, int &tsix_rna,
        vector<int> &all_xist, int &temp_sum_xr, int &temp_sum_tr, int &temp_sum_xp, int &temp_sum_tp,
        boost::random::uniform_int_distribution<> &dist, int &step, double &p8)
{
    xist_gene.push_front(xpol_ini);
    sum_xp = sum_xp - xist_gene.back() + xist_gene.front();
    xist_rna = xist_rna + (int)xist_gene.back();
    xist_gene.pop_back();
    if (p8>0) {
        collisions (dimx_xp, xist_gene, tsix_gene, sum_xp, sum_tp, dist);
    }
    tsix_gene.push_back(tpol_ini);
    sum_tp = sum_tp - tsix_gene.front() + tsix_gene.back();
    tsix_rna = tsix_rna + (int)tsix_gene.front();
    tsix_gene.pop_front();
    if (p8>0) {
        collisions (dimx_xp, xist_gene, tsix_gene, sum_xp, sum_tp, dist);
    }
    all_xist[step] = xist_rna;

    temp_sum_xr = temp_sum_xr + xist_rna;
    temp_sum_tr = temp_sum_tr + tsix_rna;
    temp_sum_xp = temp_sum_xp + sum_xp;
    temp_sum_tp = temp_sum_tp + sum_tp;
}

void mexFunction(int nhls, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

    //initiate random number generator
    struct timeval t1;
    gettimeofday(&t1, NULL);
    gen.seed(static_cast<unsigned int>(t1.tv_usec * t1.tv_sec));

	//declare variables // pointers that point at input and output variables
	mxArray *xist_pol, *tsix_pol, *xist_pol2, *tsix_pol2, *par, *const_par;
    mxArray *time_out, *xist_pol_out, *tsix_pol_out, *xist_rna_out, *tsix_rna_out,
            *xist_pol_out2, *tsix_pol_out2, *xist_rna_out2, *tsix_rna_out2, *test, *stable_switch_on_out;
	const mwSize *dims;
	double *xp, *tp, *xr, *tr, *xpr, *xp2, *tp2, *xr2, *tr2, *xpr2, *p, *cp,
            t, sil_threshold, t_diff, t_before;
    double *xpo, *tpo, *xro, *tro, *xpro, *xpo2, *tpo2, *xro2, *tro2, *xpro2,
            *time, v = 1.0/1440.0, out_step, *to, *stable_sw_on,
            k_adv_sil, k2_red_X1, k2_red_X2, k1_ind_X1, k1_ind_X2, tp_dox;
	int dimx_xp, dimx_tp, xist_rna, tsix_rna, xist_prom, xist_rna2, tsix_rna2, xist_prom2;
	int n_steps,n_steps_before, n_steps_dox, n_out_steps, n_small_steps, sec_chr, step, sel_rx;

	//associate inputs: first chromosome
    par = mxDuplicateArray(prhs[0]);
    const_par = mxDuplicateArray(prhs[1]);
	xist_pol = mxDuplicateArray(prhs[2]);
    tsix_pol = mxDuplicateArray(prhs[3]);
    xist_rna = mxGetScalar(prhs[4]);
    tsix_rna = mxGetScalar(prhs[5]);
    xist_prom = mxGetScalar(prhs[6]);

	//associate pointers
    xp = mxGetPr(xist_pol);
    tp = mxGetPr(tsix_pol);
    p = mxGetPr(par);
    cp = mxGetPr(const_par);
    t = cp[0];
    t_before =  cp[1];
    t_diff =  cp[2];
    sil_threshold = cp[3];
    out_step = cp[4];
    sec_chr = cp[5];
    k_adv_sil = cp[6];
    k2_red_X1 = cp[7];
    k2_red_X2 = cp[8];
    k1_ind_X1 = cp[9];
    k1_ind_X2 = cp[10];
    tp_dox = cp[11];

    //associate inputs: second chromosome
   if (sec_chr==1) {
        xist_pol2 = mxDuplicateArray(prhs[7]);
        tsix_pol2 = mxDuplicateArray(prhs[8]);
        xist_rna2 = mxGetScalar(prhs[9]);
        tsix_rna2 = mxGetScalar(prhs[10]);
        xist_prom2 = mxGetScalar(prhs[11]);

        xp2 = mxGetPr(xist_pol2);
        tp2 = mxGetPr(tsix_pol2);
   }

    n_out_steps = 1+ceil((t_diff+t_before)/out_step);
    n_small_steps = round(out_step/v);
    n_steps = (n_out_steps-1)*n_small_steps;
    n_steps_before = (ceil(t_before/out_step))*n_small_steps;
    n_steps_dox = (ceil(tp_dox/out_step))*n_small_steps;

	//figure out dimensions and associate output
	dims = mxGetDimensions(prhs[2]);
	dimx_xp = (int)dims[0];
    dims = mxGetDimensions(prhs[3]);
	dimx_tp = (int)dims[0];

    //output
    time_out = plhs[0] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    xist_pol_out = plhs[1] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    tsix_pol_out = plhs[2] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    xist_rna_out = plhs[3] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    tsix_rna_out = plhs[4] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    xpo = mxGetPr(xist_pol_out);
    tpo= mxGetPr(tsix_pol_out);
    xro = mxGetPr(xist_rna_out);
    tro = mxGetPr(tsix_rna_out);
    time =   mxGetPr(time_out);

    xist_pol_out2 = plhs[5] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    tsix_pol_out2 = plhs[6] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    xist_rna_out2 = plhs[7] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    tsix_rna_out2 = plhs[8] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    test = plhs[9] = mxCreateDoubleMatrix(n_out_steps, 1, mxREAL);
    stable_switch_on_out = plhs[10] = mxCreateDoubleMatrix(1, 1, mxREAL);

    xpo2 = mxGetPr(xist_pol_out2);
    tpo2= mxGetPr(tsix_pol_out2);
    xro2 = mxGetPr(xist_rna_out2);
    tro2 = mxGetPr(tsix_rna_out2);
    to = mxGetPr(test);
    stable_sw_on = mxGetPr(stable_switch_on_out);
    // Initialize
    stable_sw_on[0] = t_diff+t_before;

    //write pol mxarray in deque and sum up polymerases
    deque<int> xist_gene, tsix_gene, xist_gene2, tsix_gene2;
    int q, sum_xp=0, sum_tp=0, sum_xp2 = 0, sum_tp2 = 0;

    for (q=0; q<dimx_xp; q++)
        {
            xist_gene.push_back (xp[q]);
            sum_xp = sum_xp + xp[q];
             if (sec_chr==1) {
                xist_gene2.push_back (xp2[q]);
                sum_xp2 = sum_xp2 + xp2[q];
             }
        }
    for (q=0; q<dimx_tp; q++)
        {
            tsix_gene.push_back (tp[q]);
            sum_tp = sum_tp + tp[q];
            if (sec_chr==1) {
                tsix_gene2.push_back (tp2[q]);
                 sum_tp2 = sum_tp2 + tp2[q];
            }
        }

    bool xpol_ini=0, tpol_ini=0, xpol_ini2=0, tpol_ini2=0;
    boost::random::uniform_int_distribution<> dist(0, 1);
    boost::random::uniform_real_distribution<> dist_real(0.0, 1.0);
    double t_next, r, sum_rx;
    vector<double> rx(16), cum_rx(16);
    vector<int> all_xist(n_steps), all_xist2(n_steps);

    int sil_xa = 0, sil_xa2 = 0, xa_sil_steps = round(p[11]/v),
            sil_tsix = 0, sil_tsix2 = 0,
	    sil_state_xa = 0, sil_state_xa2 = 0, sil_state_tsix = 0, sil_state_tsix2 = 0, stst = 0;
    double xa = 1, xa2 = 0;
    if (sec_chr==1) {xa2 = 1;};
    double temp_sum_xa = xa + xa2;
    int m = 0, b = 0, temp_sum_xr = 0, temp_sum_tr = 0, temp_sum_xp = 0, temp_sum_tp = 0,
            temp_sum_xr2 = 0, temp_sum_tr2 = 0, temp_sum_xp2 = 0, temp_sum_tp2 = 0;
    all_xist[0] = xist_rna;
    xro[b] = xist_rna;
    tro[b] = tsix_rna;
    xpo[b] = sum_xp;
    tpo[b] = sum_tp;
    to[b] = temp_sum_xa;
    if (sec_chr==1) {
        all_xist2[0] = xist_rna2;
        xro2[b] = xist_rna2;
        tro2[b] = tsix_rna2;
        xpo2[b] = sum_xp2;
        tpo2[b] = sum_tp2;
    }
    //while (act_time<t_max)
    for (step=0; step<n_steps; step++)
        {
			// Identify first timepoint at which stable state with 1 Xa and 1 Xi is reached
			if (xist_rna>=sil_threshold & xa ==0 & sil_tsix ==1 & xist_rna2<sil_threshold & xa2==1 & sil_tsix2==0 & stst==0) {stable_sw_on[0] = step*v; stst=1;}
			else { if (xist_rna2>=sil_threshold & xa2 ==0 & sil_tsix2 ==1 & xist_rna<sil_threshold & xa==1 & sil_tsix==0 & stst==0) {stable_sw_on[0] = step*v; stst=1;}
			}
			//Silencing
			//If Tsix silencing has reached final state, Tsix is silenced //p[12]=parameter that determines # of intermediate states until silencing of tsix occurs
			if (sil_state_tsix>=p[12]&xist_rna>=sil_threshold) {sil_tsix=1;}
			//If silencing has not reached the final state but Xist rna level is below threshold system falls back into initial tsix silencing state
			else {
			if (xist_rna<sil_threshold&sil_state_tsix<p[12]) {sil_state_tsix=0; sil_tsix=0;}
			}
			//If XA silencing has reached final state, XA is silenced //p[11]=parameter that determines # of intermediate states until silencing of xa occurs
			if (sil_state_xa>=p[11]&xist_rna>=sil_threshold) {xa=0;}
			//If silencing has not reached the final state but Xist rna level is below threshold system falls back into initial xa silencing state
			else {
			if (xist_rna<sil_threshold&sil_state_xa<p[11]) {sil_state_xa=0; xa=1;}
			}
			if (p[6]==0&tsix_gene.front()==1)
		        	{
		            	xist_prom = 2;
				}
		        if (sec_chr==1) {
		     	if (sil_state_tsix2>=p[12]&xist_rna2>=sil_threshold) {sil_tsix2=1;}
			else {
			if (xist_rna2<sil_threshold&sil_state_tsix2<p[12]) {sil_state_tsix2=0; sil_tsix2=0;}
			}
			if (sil_state_xa2>=p[11]&xist_rna2>=sil_threshold) {xa2=0;}
			else {
			if (xist_rna2<sil_threshold&sil_state_xa2<p[11]) {sil_state_xa2=0; xa2=1;}
			}
			if (p[6]==0&tsix_gene2.front()==1)
				{
		                	xist_prom2 = 2;
				}
		        }

			//Gillespie reactions
	        t_next = 0;

	        xpol_ini=0, tpol_ini=0, xpol_ini2=0, tpol_ini2=0;
	        while (t_next<v) {
				//Xist initiation
	            //before induction of diff k1 10 times reduced
	            if (step<n_steps_before & step<n_steps_dox) {rx[0] = p[0]*0*(xa+xa2)*(xist_prom==1)*(tsix_gene.front()==0);}
	            else {
					if (step<n_steps_dox) {rx[0] = p[0]*(xa+xa2)*(xist_prom==1)*(tsix_gene.front()==0);}
					// if Xist is induced with Dox (k1_ind_X1>0) the Xist initiation rate becomes independent of the XA concentration
					else {
						if (k1_ind_X1==0) {rx[0] = (p[0]*(step>n_steps_before))*(xa+xa2)*(xist_prom==1)*(tsix_gene.front()==0);}
						// If X1 is induced with dox, the Xist initiation rate becomes independent of xa+xa2
						else {rx[0] = (p[0]*(step>n_steps_before)+p[0]*k1_ind_X1)*(xist_prom==1)*(tsix_gene.front()==0);}
					}
				}
	            rx[1] = (p[1]-k2_red_X1*p[1])*(sil_tsix<1); //Tsix initiation
	            rx[2] = p[7]*(xist_prom==2)*(tsix_gene.front()==0); // turning on the Xist promoter
	            rx[3] = p[3]*xist_rna; //Xist degradation
	            rx[4] = p[4]*tsix_rna; //Tsix degradation

	            //second chr
	            //Xist initiation
	            if (step<n_steps_before & step<n_steps_dox) {rx[5] = (sec_chr==1)*p[0]*0*(xa+xa2)*(xist_prom2==1)*(tsix_gene2.front()==0);}
	            else {
					if (step<n_steps_dox) {rx[5] = (sec_chr==1)*p[0]*(xa+xa2)*(xist_prom2==1)*(tsix_gene2.front()==0);}
					else {
						if (k1_ind_X2==0) {rx[5] = (sec_chr==1)*(p[0]*(step>n_steps_before))*(xa+xa2)*(xist_prom2==1)*(tsix_gene2.front()==0);}
						// If X2 is induced with dox, the Xist initiation rate becomes independent of xa+xa2
						else {rx[5] = (sec_chr==1)*(p[0]*(step>n_steps_before)+p[0]*k1_ind_X2)*(xist_prom2==1)*(tsix_gene2.front()==0);}
					}
				}
	            rx[6] = (sec_chr==1)*(p[1]-k2_red_X2*p[1])*(sil_tsix2<1);   //Tsix initiation
	            rx[7] = (sec_chr==1)*(p[7]*(xist_prom2==2))*(tsix_gene2.front()==0); // turning on the Xist promoter
	            rx[8] = (sec_chr==1)*p[3]*xist_rna2; //Xist degradation
	            rx[9] = (sec_chr==1)*p[4]*tsix_rna2; //Tsix degradation

		    //transition between silencing states as Gillespie reaction
		    ////chromosome 1: now k_adv_sil= rate for transition between silencing states
		    rx[10] = k_adv_sil*(xist_rna>=sil_threshold)*(sil_state_xa<p[11]|sil_state_tsix<p[12]); //XA and Tsix silencing advance one state
		    ////if system is in last silencing state of the xa but xist rna fall below the threshold the silencing state can go back to the initial state with rate k[13]
		    rx[11] = p[13]*(xist_rna<sil_threshold)*(sil_state_xa>=p[11]);
		    //if system is in last silencing state of tsix but xist rna fall below the threshold the silencing state can go back to the initial state with rate k[14]
		    rx[12] = p[14]*(xist_rna<sil_threshold)*(sil_state_tsix>=p[12]);
		    ////second chromosome
		    rx[13] = (sec_chr==1)*k_adv_sil*(xist_rna2>=sil_threshold)*(sil_state_xa2<p[11]|sil_state_tsix2<p[12]);
		    rx[14] = (sec_chr==1)*p[13]*(xist_rna2<sil_threshold)*(sil_state_xa2>=p[11]);
		    rx[15] = (sec_chr==1)*p[14]*(xist_rna2<sil_threshold)*(sil_state_tsix2>=p[12]);

            sum_rx = 0, sel_rx = 100;
            for (q=0; q<rx.size(); q++) {
                sum_rx = sum_rx + rx[q];
                cum_rx[q] = sum_rx;
            }

            if (sum_rx>0) {
                t_next = t_next + (-log(dist_real(gen)))/(sum_rx);
            } else {t_next=1;}

            if (t_next<v){
                // decide which reaction occurs
                r = dist_real(gen);
                for (q=0; q<rx.size(); q++) {
                    if (r<(cum_rx[q]/sum_rx)){
                        sel_rx = q;
                        break;
                    }
                }
                //execute reaction

                switch (sel_rx){
                    case 0:
                        xpol_ini = 1; break;
                    case 1:
                        tpol_ini = 1; break;
                    case 2:
                        xist_prom = 1; break;
                    case 3:
                        xist_rna--;
                         break;
                    case 4:
                        tsix_rna--; break;

                    //second chr

                   case 5:
                        xpol_ini2 = 1; break;
                    case 6:
                        tpol_ini2 = 1; break;
                    case 7:
                        xist_prom2 = 1; break;
                    case 8:
                        xist_rna2--;
                        break;
                    case 9:
                        tsix_rna2--; break;
					case 10:
						if (sil_state_xa<p[11]) {sil_state_xa++;}
						if (sil_state_tsix<p[12]) {sil_state_tsix++;}
						break;
					case 11:
						sil_state_xa=0; xa=1; break;
					case 12:
						sil_state_tsix=0; sil_tsix=0; break;
					case 13:
						if (sil_state_xa2<p[11]) {sil_state_xa2++;}
						if (sil_state_tsix2<p[12]) {sil_state_tsix2++;}
						break;
					case 14:
						sil_state_xa2=0; xa2=1; break;
					case 15:
						sil_state_tsix2=0; sil_tsix2=0; break;
					default:
                        break;
                }
            }
        }
        temp_sum_xa = temp_sum_xa + xa + xa2;
        //elongation and collisions
        elongation (dimx_xp, xist_gene, tsix_gene, sum_xp, sum_tp,
            xpol_ini, tpol_ini, xist_rna, tsix_rna,
            all_xist, temp_sum_xr, temp_sum_tr, temp_sum_xp, temp_sum_tp, dist, step, p[8]);

        //second chromosome
        if (sec_chr==1) {
            elongation (dimx_xp, xist_gene2, tsix_gene2, sum_xp2, sum_tp2,
            xpol_ini2, tpol_ini2, xist_rna2, tsix_rna2,
            all_xist2, temp_sum_xr2, temp_sum_tr2, temp_sum_xp2, temp_sum_tp2, dist, step, p[8]);

        }

        //write out smoothed variable
        m++;
        if (m==n_small_steps) {
            b++;
            xro[b] = temp_sum_xr/n_small_steps;
            tro[b] = temp_sum_tr/n_small_steps;
            xpo[b] = temp_sum_xp/n_small_steps;
            tpo[b] = temp_sum_tp/n_small_steps;
            to[b] = temp_sum_xa/n_small_steps;
            time[b] = step*v;
            temp_sum_xr = 0;
            temp_sum_tr = 0;
            temp_sum_xp = 0;
            temp_sum_tp = 0;
            temp_sum_xa = 0;
            if (sec_chr==1) {
                xro2[b] = temp_sum_xr2/n_small_steps;
                tro2[b] = temp_sum_tr2/n_small_steps;
                xpo2[b] = temp_sum_xp2/n_small_steps;
                tpo2[b] = temp_sum_tp2/n_small_steps;
                temp_sum_xr2 = 0;
                temp_sum_tr2 = 0;
                temp_sum_xp2 = 0;
                temp_sum_tp2 = 0;
            }
            m = 0;
        }

        }
    if (m>0) {
            b++;
            xro[b] = temp_sum_xr/m;
            tro[b] = temp_sum_tr/m;
            xpo[b] = temp_sum_xp/m;
            tpo[b] = temp_sum_tp/m;
            to[b] = temp_sum_xa/n_small_steps;
            time[b] = step*v;
            if (sec_chr==1) {
                xro2[b] = temp_sum_xr2/m;
                tro2[b] = temp_sum_tr2/m;
                xpo2[b] = temp_sum_xp2/m;
                tpo2[b] = temp_sum_tp2/m;
            }
    }
   	}


text:
/*
 * Copyright (c) 2016-2018 Nicholas Corgan (n.corgan@gmail.com)
 *
 * Distributed under the MIT License (MIT) (See accompanying file LICENSE.txt
 * or copy at http://opensource.org/licenses/MIT)
 */

#include "private_exports.hpp"

#include "database/enum_conversions.hpp"
#include "database/database_common.hpp"

#include "types/rng.hpp"

#include <pkmn/config.hpp>
#include <pkmn/utils/paths.hpp>

#include <boost/filesystem.hpp>

#ifdef PKMN_ENABLE_QT
#include <QImage>
#include <QString>
#endif

#ifdef PKMN_PLATFORM_WIN32
#    include <windows.h>
#    define GET_PROCESS_ID() GetCurrentProcessId()
#else
#    include <unistd.h>
#    define GET_PROCESS_ID() getpid()
#endif

namespace fs = boost::filesystem;

namespace pkmn { namespace priv {

    int game_name_to_generation(
        const std::string& game
    )
    {
        return pkmn::database::game_name_to_generation(game);
    }

    int game_enum_to_generation(pkmn::e_game game)
    {
        return pkmn::database::game_enum_to_generation(game);
    }

    void initialize_database_connection()
    {
        pkmn::database::initialize_connection();
    }

    std::string convert_image_for_mono(
        const std::string& filepath
    )
    {
#ifdef PKMN_ENABLE_QT
        static const fs::path TMP_DIR(pkmn::get_tmp_dir());

        QImage input_image(QString::fromStdString(filepath));
        QImage output_image = input_image.convertToFormat(QImage::Format_ARGB32);

        fs::path output_path(TMP_DIR);

        std::string output_filename(std::to_string(GET_PROCESS_ID()));
        output_filename.append(std::to_string(pkmn::rng<uint32_t>().rand()));
        output_filename.append(".png");

        output_path /= output_filename;

        output_image.save(QString::fromStdString(output_path.string()));

        return output_path.string();
#else
        return filepath;
#endif
    }

    pkmn::e_item ball_to_item(pkmn::e_ball ball)
    {
        return pkmn::database::ball_to_item(ball);
    }

    pkmn::e_ball item_to_ball(pkmn::e_item item)
    {
        return pkmn::database::item_to_ball(item);
    }
}}

ext: cpp

max_stars_repo_name: ncorgan/libpkmn

max_issues_repo_name: PMArkive/libpkmn

max_forks_repo_name: PMArkive/libpkmn


text:
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/module.hpp>
#include <boost/python/wrapper.hpp>

#include <rtamt_stl_library/stl_node.h>

using namespace boost::python;
using namespace stl_library;

struct StlNodeWrap : StlNode, wrapper<StlNode>
{
    Sample update() {
        return this->get_override("update")();
    }

    void addNewInput(int i, Sample msg) {
        this->get_override("addNewInput")(i, msg);
    }
};

BOOST_PYTHON_MODULE(stl_node)
{
    class_<StlNodeWrap, boost::noncopyable>("StlNode")
        .def("update", pure_virtual(&StlNode::update))
    ;
}

max_stars_repo_name: BentleyJOakes/rtamt

max_issues_repo_name: BentleyJOakes/rtamt

max_forks_repo_name: BentleyJOakes/rtamt


text:
//==================================================================================================
/*!

  Copyright 2016 NumScale SAS

  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
*/
//==================================================================================================
#include <boost/simd/function/scalar/is_imag.hpp>
#include <scalar_test.hpp>
#include <boost/simd/detail/dispatch/meta/as_integer.hpp>
#include <boost/simd/constant/inf.hpp>
#include <boost/simd/constant/minf.hpp>
#include <boost/simd/constant/mone.hpp>
#include <boost/simd/constant/nan.hpp>
#include <boost/simd/constant/one.hpp>
#include <boost/simd/constant/zero.hpp>
#include <boost/simd/constant/quarter.hpp>
#include <boost/simd/constant/half.hpp>
#include <boost/simd/constant/two.hpp>
#include <boost/simd/logical.hpp>

STF_CASE_TPL (" is_imag real",  STF_IEEE_TYPES)
{
  namespace bs = boost::simd;
  namespace bd = boost::dispatch;
  using bs::is_imag;
  using r_t = decltype(is_imag(T()));

  // return type conformity test
  STF_TYPE_IS(r_t, bs::logical<T>);

  // specific values tests
#ifndef BOOST_SIMD_NO_INVALIDS
  STF_EQUAL(is_imag(bs::Inf<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Minf<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Nan<T>()), r_t(false));
#endif
  STF_EQUAL(is_imag(-bs::Zero<T>()), r_t(true));
  STF_EQUAL(is_imag(bs::Half<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Mone<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::One<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Quarter<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Two<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Zero<T>()), r_t(true));
}

STF_CASE_TPL (" is_imag signed_int",  STF_SIGNED_INTEGRAL_TYPES)
{
  namespace bs = boost::simd;
  namespace bd = boost::dispatch;
  using bs::is_imag;
   using r_t = decltype(is_imag(T()));

  // return type conformity test
  STF_TYPE_IS(r_t, bs::logical<T>);

  // specific values tests
  STF_EQUAL(is_imag(bs::Mone<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::One<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Two<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Zero<T>()), r_t(true));
}

STF_CASE_TPL (" is_imag unsigned_int",  STF_UNSIGNED_INTEGRAL_TYPES)
{
  namespace bs = boost::simd;
  namespace bd = boost::dispatch;
  using bs::is_imag;
  using r_t = decltype(is_imag(T()));

  // return type conformity test
  STF_TYPE_IS(r_t, bs::logical<T>);

  // specific values tests
  STF_EQUAL(is_imag(bs::One<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Two<T>()), r_t(false));
  STF_EQUAL(is_imag(bs::Zero<T>()), r_t(true));
}

STF_CASE ("is_imag bool")
{
  namespace bs = boost::simd;
  namespace bd = boost::dispatch;
  using bs::is_imag;
   using r_t = decltype(is_imag(bool()));

  // return type conformity test
   STF_TYPE_IS(r_t, bool);

   // specific values tests
  STF_EQUAL(is_imag(true), false);
  STF_EQUAL(is_imag(false), true);
}


text:
///////////////////////////////////////////////////////////////////////////////
// ars::include.hpp                                                          //
//                                                                           //
//  Copyright 2009 Erwann Rogard. Distributed under the Boost                //
//  Software License, Version 1.0. (See accompanying file                    //
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)         //
///////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_STATISTICS_DETAIL_ARS_INCLUDE_HPP_ER_2009
#define BOOST_STATISTICS_DETAIL_ARS_INCLUDE_HPP_ER_2009

#include <boost/ars/constant.hpp>
#include <boost/ars/proposal_sampler.hpp>
#include <boost/ars/sampler.hpp>

#endif

ext: hpp


text:
// Copyright (C) 2006  Davis E. King (davis@dlib.net)
// License: Boost Software License   See LICENSE.txt for the full license.

#include <sstream>
#include <string>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <dlib/rand.h>
#include <dlib/compress_stream.h>
#include <dlib/hash.h>

#include "tester.h"

namespace
{

    using namespace test;
    using namespace dlib;
    using namespace std;

    logger dlog("test.rand");

    void check_bpp (
        const std::string str
    )
    {
        istringstream rdata;
        ostringstream sout;
        rdata.str(str);
        double compressed_size;
        compress_stream::kernel_1a cs1;
        compress_stream::kernel_2a cs2;

        compress_stream_kernel_1<
            entropy_encoder_model_kernel_5<257,entropy_encoder::kernel_1a,4000000,4>,
            entropy_decoder_model_kernel_5<257,entropy_decoder::kernel_1a,4000000,4>,
            crc32::kernel_1a
            > cs3;

        print_spinner();

        rdata.clear();
        rdata.seekg(0);
        sout.clear();
        sout.str("");
        cs1.compress(rdata,sout);
        compressed_size = sout.str().size();
        compressed_size *= 8;
        compressed_size /= str.size();
        DLIB_TEST_MSG(compressed_size >= 8, "order 0 bps: " << compressed_size);
        dlog << LINFO << "order 0: " << compressed_size;

        print_spinner();

        rdata.clear();
        rdata.seekg(0);
        sout.clear();
        sout.str("");
        cs2.compress(rdata,sout);
        compressed_size = sout.str().size();
        compressed_size *= 8;
        compressed_size /= str.size();
        DLIB_TEST_MSG(compressed_size >= 8, "order 1 bps: " << compressed_size);
        dlog << LINFO << "order 1: " << compressed_size;

        print_spinner();

        rdata.clear();
        rdata.seekg(0);
        sout.clear();
        sout.str("");
        cs3.compress(rdata,sout);
        compressed_size = sout.str().size();
        compressed_size *= 8;
        compressed_size /= str.size();
        DLIB_TEST_MSG(compressed_size >= 8, "order 4 bps: " << compressed_size);
        dlog << LINFO << "order 4: " << compressed_size;

    }

    template <
        typename rand
        >
    void rand_test (
    )
    /*!
        requires
            - rand is an implementation of rand/rand_kernel_abstract.h
              is instantiated with int
        ensures
            - runs tests on rand for compliance with the specs
    !*/
    {

        ostringstream seed;
        seed << (unsigned int)time(0);

        ostringstream sout;

        rand r, r2;
        DLIB_TEST(r.get_seed() == "");
        r.set_seed(seed.str());

        DLIB_TEST(r.get_seed() == seed.str());
        r.clear();
        DLIB_TEST(r.get_seed() == "");
        swap(r,r2);
        DLIB_TEST(r.get_seed() == "");
        r.set_seed(seed.str());
        DLIB_TEST(r.get_seed() == seed.str());
        swap(r,r2);
        DLIB_TEST(r2.get_seed() == seed.str());
        DLIB_TEST(r.get_seed() == "");
        swap(r,r2);
        DLIB_TEST(r.get_seed() == seed.str());
        DLIB_TEST(r2.get_seed() == "");

        print_spinner();
        unsigned long size = 100000;
        for (unsigned long i = 0; i < size; ++i)
        {
            uint32 ch = r.get_random_32bit_number();
            sout.write((char*)&ch,4);
        }

        check_bpp(sout.str());
        sout.clear();
        sout.str("");

        print_spinner();
        for (unsigned long i = 0; i < size; ++i)
        {
            uint16 ch = r.get_random_16bit_number();
            sout.write((char*)&ch,2);
        }

        check_bpp(sout.str());
        sout.clear();
        sout.str("");

        print_spinner();
        for (unsigned long i = 0; i < size; ++i)
        {
            unsigned char ch = r.get_random_8bit_number();
            sout.write((char*)&ch,1);
        }

        check_bpp(sout.str());
        sout.clear();
        sout.str("");

        // make sure the things can serialize right
        {
            r.clear();
            r2.clear();

            for (int i =0; i < 1000; ++i)
            {
                r.get_random_32bit_number();
                r.get_random_gaussian();
            }

            ostringstream sout;
            serialize(r, sout);

            istringstream sin(sout.str());
            deserialize(r2, sin);

            for (int i =0; i < 1000; ++i)
            {
                DLIB_TEST(r.get_random_32bit_number() == r2.get_random_32bit_number());
                DLIB_TEST(std::abs(r.get_random_gaussian() - r2.get_random_gaussian()) < 1e-14);
            }
        }

        // make sure calling clear() and set_seed("") do the same thing
        {
            r.clear();
            r2.set_seed("");
            rand r3;

            DLIB_TEST(r.get_seed() == r2.get_seed());
            DLIB_TEST(r.get_seed() == r3.get_seed());

            for (int i =0; i < 1000; ++i)
            {
                const uint32 num1 = r.get_random_32bit_number();
                const uint32 num2 = r2.get_random_32bit_number();
                const uint32 num3 = r3.get_random_32bit_number();
                DLIB_TEST( num1 == num2);
                DLIB_TEST( num1 == num3);
            }
        }

    }

    template <typename rand_type>
    void test_normal_numbers(
        rand_type& rnd
    )
    {
        print_spinner();
        dlog << LINFO << "test normality";
        double cnt1 = 0; // num <= -1.2
        double cnt2 = 0; // num <= -0.5
        double cnt3 = 0; // num <= 0
        double cnt4 = 0; // num <= 0.5
        double cnt5 = 0; // num <= 1.2

        const unsigned long total = 1000000;
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = rnd.get_random_gaussian();
            if (r <= -1.2) cnt1 += 1;
            if (r <= -0.5) cnt2 += 1;
            if (r <=  0)   cnt3 += 1;
            if (r <=  0.5) cnt4 += 1;
            if (r <=  1.2) cnt5 += 1;
        }

        cnt1 /= total;
        cnt2 /= total;
        cnt3 /= total;
        cnt4 /= total;
        cnt5 /= total;

        dlog << LINFO << "cnt1: "<< cnt1;
        dlog << LINFO << "cnt2: "<< cnt2;
        dlog << LINFO << "cnt3: "<< cnt3;
        dlog << LINFO << "cnt4: "<< cnt4;
        dlog << LINFO << "cnt5: "<< cnt5;

        DLIB_TEST(std::abs(cnt1 - 0.11507) < 0.001);
        DLIB_TEST(std::abs(cnt2 - 0.30854) < 0.001);
        DLIB_TEST(std::abs(cnt3 - 0.5)     < 0.001);
        DLIB_TEST(std::abs(cnt4 - 0.69146) < 0.001);
        DLIB_TEST(std::abs(cnt5 - 0.88493) < 0.001);

    }

    void test_gaussian_random_hash()
    {
        print_spinner();
        dlog << LINFO << "test_gaussian_random_hash()";
        double cnt1 = 0; // num <= -1.2
        double cnt2 = 0; // num <= -0.5
        double cnt3 = 0; // num <= 0
        double cnt4 = 0; // num <= 0.5
        double cnt5 = 0; // num <= 1.2

        const unsigned long total = 1000000;
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = gaussian_random_hash(i,0,0);
            if (r <= -1.2) cnt1 += 1;
            if (r <= -0.5) cnt2 += 1;
            if (r <=  0)   cnt3 += 1;
            if (r <=  0.5) cnt4 += 1;
            if (r <=  1.2) cnt5 += 1;
        }
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = gaussian_random_hash(0,i,0);
            if (r <= -1.2) cnt1 += 1;
            if (r <= -0.5) cnt2 += 1;
            if (r <=  0)   cnt3 += 1;
            if (r <=  0.5) cnt4 += 1;
            if (r <=  1.2) cnt5 += 1;
        }
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = gaussian_random_hash(0,0,i);
            if (r <= -1.2) cnt1 += 1;
            if (r <= -0.5) cnt2 += 1;
            if (r <=  0)   cnt3 += 1;
            if (r <=  0.5) cnt4 += 1;
            if (r <=  1.2) cnt5 += 1;
        }

        cnt1 /= total*3;
        cnt2 /= total*3;
        cnt3 /= total*3;
        cnt4 /= total*3;
        cnt5 /= total*3;

        dlog << LINFO << "cnt1: "<< cnt1;
        dlog << LINFO << "cnt2: "<< cnt2;
        dlog << LINFO << "cnt3: "<< cnt3;
        dlog << LINFO << "cnt4: "<< cnt4;
        dlog << LINFO << "cnt5: "<< cnt5;

        DLIB_TEST(std::abs(cnt1 - 0.11507) < 0.001);
        DLIB_TEST(std::abs(cnt2 - 0.30854) < 0.001);
        DLIB_TEST(std::abs(cnt3 - 0.5)     < 0.001);
        DLIB_TEST(std::abs(cnt4 - 0.69146) < 0.001);
        DLIB_TEST(std::abs(cnt5 - 0.88493) < 0.001);
    }

    void test_uniform_random_hash()
    {
        print_spinner();
        dlog << LINFO << "test_uniform_random_hash()";
        double cnt1 = 0; // num <= 0.2
        double cnt2 = 0; // num <= 0.4
        double cnt3 = 0; // num <= 0.6
        double cnt4 = 0; // num <= 0.8
        double cnt5 = 0; // num <= 1.0

        double min_val = 10;
        double max_val = 0;

        const unsigned long total = 1000000;
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = uniform_random_hash(i,0,0);
            min_val = min(r,min_val);
            max_val = max(r,max_val);

            if (r <=  0.2) cnt1 += 1;
            if (r <=  0.4) cnt2 += 1;
            if (r <=  0.6) cnt3 += 1;
            if (r <=  0.8) cnt4 += 1;
            if (r <=  1.0) cnt5 += 1;
        }
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = uniform_random_hash(0,i,0);
            min_val = min(r,min_val);
            max_val = max(r,max_val);

            if (r <=  0.2) cnt1 += 1;
            if (r <=  0.4) cnt2 += 1;
            if (r <=  0.6) cnt3 += 1;
            if (r <=  0.8) cnt4 += 1;
            if (r <=  1.0) cnt5 += 1;
        }
        for (unsigned long i = 0; i < total; ++i)
        {
            const double r = uniform_random_hash(0,0,i);
            min_val = min(r,min_val);
            max_val = max(r,max_val);

            if (r <=  0.2) cnt1 += 1;
            if (r <=  0.4) cnt2 += 1;
            if (r <=  0.6) cnt3 += 1;
            if (r <=  0.8) cnt4 += 1;
            if (r <=  1.0) cnt5 += 1;
        }

        cnt1 /= total*3;
        cnt2 /= total*3;
        cnt3 /= total*3;
        cnt4 /= total*3;
        cnt5 /= total*3;

        dlog << LINFO << "cnt1: "<< cnt1;
        dlog << LINFO << "cnt2: "<< cnt2;
        dlog << LINFO << "cnt3: "<< cnt3;
        dlog << LINFO << "cnt4: "<< cnt4;
        dlog << LINFO << "cnt5: "<< cnt5;
        dlog << LINFO << "min_val: "<< min_val;
        dlog << LINFO << "max_val: "<< max_val;

        DLIB_TEST(std::abs(cnt1 - 0.2) < 0.001);
        DLIB_TEST(std::abs(cnt2 - 0.4) < 0.001);
        DLIB_TEST(std::abs(cnt3 - 0.6) < 0.001);
        DLIB_TEST(std::abs(cnt4 - 0.8) < 0.001);
        DLIB_TEST(std::abs(cnt5 - 1.0) < 0.001);
        DLIB_TEST(std::abs(min_val - 0.0) < 0.001);
        DLIB_TEST(std::abs(max_val - 1.0) < 0.001);
    }

    void test_get_integer()
    {

        print_spinner();
        dlib::rand rnd;

        int big_val = 0;
        int small_val = 0;

        const long long maxval = (((unsigned long long)1)<<62) + (((unsigned long long)1)<<61);
        for (int i = 0; i < 10000000; ++i)
        {
            if (rnd.get_integer(maxval) > maxval/2)
                ++big_val;
            else
                ++small_val;
        }

        // make sure there isn't any funny bias
        DLIB_TEST(std::abs(big_val/(double)small_val - 1) < 0.001);

        //cout << big_val/(double)small_val << endl;

    }

    class rand_tester : public tester
    {
    public:
        rand_tester (
        ) :
            tester ("test_rand",
                    "Runs tests on the rand component.")
        {}

        void perform_test (
        )
        {
            dlog << LINFO << "testing kernel_1a";
            rand_test<dlib::rand>();
            rand_test<dlib::rand>();

            dlib::rand rnd;
            test_normal_numbers(rnd);
            test_gaussian_random_hash();
            test_uniform_random_hash();
            test_get_integer();
        }
    } a;

}


text:
/* Copyright (c) 2010-2017 Stanford University
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR(S) DISCLAIM ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL AUTHORS BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

#include <fcntl.h>
#include <stdarg.h>
#include <execinfo.h>
#include <signal.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdexcept>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Weffc++"
#include <boost/lexical_cast.hpp>
#pragma GCC diagnostic pop

#include "Cycles.h"
#include "LogCabinLogger.h"
#include "Logger.h"
#include "ShortMacros.h"
#include "ThreadId.h"
#include "TimeTrace.h"
#include "Unlock.h"
#include "Util.h"

namespace RAMCloud {

/**
 * Friendly names for each #LogLevel value.
 * Keep this in sync with the LogLevel enum.
 */
static const char* logLevelNames[] = {"(none)", "ERROR", "WARNING",
                                      "NOTICE", "DEBUG"};

static_assert(unsafeArrayLength(logLevelNames) == NUM_LOG_LEVELS,
              "logLevelNames size does not match NUM_LOG_LEVELS");
/**
 * Friendly names for each #LogModule value.
 * Keep this in sync with the LogModule enum.
 */
static const char* logModuleNames[] = {
    "default",
    "transport",
    "externalStorage"
};

static_assert(unsafeArrayLength(logModuleNames) == NUM_LOG_MODULES,
              "logModuleNames size does not match NUM_LOG_MODULES");

/**
 * Create a new debug logger; messages will go to stderr by default. Should
 * not be called outside this class except during unit testing.
 * \param[in] level
 *      Messages for all modules at least as important as \a level will be
 *      logged.
 */
Logger::Logger(LogLevel level)
    : fd(2)
    , mustCloseFd(false)
    , mutex("Logger::mutex")
    , collapseMap()
    , collapseIntervalMs(DEFAULT_COLLAPSE_INTERVAL)
    , maxCollapseMapSize(DEFAULT_COLLAPSE_MAP_LIMIT)
    , nextCleanTime({0, 0})
    , logDataAvailable()
    , bufferSpaceFreed()
    , bufferSize(1000000)
    , messageBuffer(new char[bufferSize])
    , nextToInsert(0)
    , nextToPrint(0)
    , discardedEntries(0)
    , printThread()
    , printThreadExit(false)
    , testingBufferSize(0)
    , testingNoNotify(false)
    , testingLogTime(NULL)
    , testingStallPrintThread(false)
{
    setLogLevels(level);

    // Touch every page in the buffer to make sure that memory has been
    // allocated (don't want to take page faults during log operations).
    for (int i = 0; i < bufferSize; i += 1000) {
        messageBuffer[i] = 'x';
    }

    // Start the print thread.
    printThread.construct(printThreadMain, this);

#if ENABLE_LOGCABIN
    LogCabinLogger::setup(level);
#endif
}

/**
 * Destructor for debug logs.
 */
Logger::~Logger()
{
    // Exit the print thread.
    {
        Lock lock(mutex);
        printThreadExit = true;
        logDataAvailable.notify_one();
    }
    printThread->join();

    // No lock needed: the print thread is finished
    if (mustCloseFd)
        close(fd);
    delete[] messageBuffer;
}

/**
 * Return the singleton shared instance that is normally used for logging.
 */
Logger&
Logger::get()
{
    // Use static local variable to achieve efficient thread-safe lazy
    // initialization. If multiple threads attempt to initialize sharedLogger
    // concurrently, the initialization is guaranteed to occur exactly once.
    static Logger sharedLogger;
    return sharedLogger;
}

/**
 * Arrange for future log messages to go to a particular file.
 *
 * This should only be invoked when one is sure that printThread is not
 * running (e.g. before any logMessage() invocations) due to a file
 * descriptor race that can occur (RAM-839).
 * \param path
 *      The pathname for the log file, which may or may not exist already.
 * \param truncate
 *      True means the log file should be truncated if it already exists;
 *      false means an existing log file is retained, and new messages are
 *      appended.
 */
void
Logger::setLogFile(const char* path, bool truncate)
{
    Lock lock(mutex);
    int newFd = open(path, O_CREAT | O_WRONLY | (truncate ? O_TRUNC : 0),
            0666);
    if (newFd < 0) {
        throw Exception(HERE,
                        format("couldn't open log file '%s'", path),
                        errno);
    }
    lseek(newFd, 0, SEEK_END);
    if (mustCloseFd)
        close(fd);
    fd = newFd;
    mustCloseFd = true;
}

/**
 * Arrange for future log messages to go to a particular file descriptor.
 * \param newFd
 *      The file descriptor for the log file, which must already have
 *      been opened by the caller. This file descriptor belongs to the caller,
 *      meaning that we will never close it.
 */
void
Logger::setLogFile(int newFd)
{
    Lock lock(mutex);
    if (mustCloseFd)
        close(fd);
    fd = newFd;
    mustCloseFd = false;
}

/**
 * Set the log level for a particular module.
 * \param[in] module
 *      The module whose level to set.
 * \param[in] level
 *      Messages for \a module at least as important as \a level will be
 *      logged.
 */
void
Logger::setLogLevel(LogModule module, LogLevel level)
{
    Lock lock(mutex);
    logLevels[module] = level;
}

/**
 * Set the log level for a particular module.
 * \param[in] module
 *      The module whose level to set.
 * \param[in] level
 *      Messages for \a module at least as important as \a level will be
 *      logged. This will be clamped to a valid LogLevel if it is out of range.
 */
void
Logger::setLogLevel(LogModule module, int level)
{
    // No lock needed: doesn't access Logger object.
    if (level < 0)
        level = 0;
    else if (level >= NUM_LOG_LEVELS)
        level = NUM_LOG_LEVELS - 1;
    setLogLevel(module, static_cast<LogLevel>(level));
}

/**
 * Set the log level for a particular module.
 * \param[in] module
 *      The module whose level to set.  Given as a string from #logModuleNames.
 * \param[in] level
 *      Messages for \a module at least as important as \a level will be
 *      logged. This will be clamped to a valid LogLevel if it is out of range.
 *      Given as a string from #logLevelNames or as a decimal number indicating
 *      that level's index in the array.
 */
void
Logger::setLogLevel(string module, string level)
{
    // No lock needed: doesn't access Logger object.
    int moduleIndex = 0;
    for (; moduleIndex < NUM_LOG_MODULES; ++moduleIndex) {
        if (module == logModuleNames[moduleIndex])
            break;
    }
    if (moduleIndex == NUM_LOG_MODULES) {
        LOG(WARNING, "Ignoring bad log module name: %s", module.c_str());
        return;
    }
    int moduleLevel;
// The uninitialized variable is inside boost::lexical_cast.
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
    try {
        moduleLevel = boost::lexical_cast<int>(level);
    } catch (boost::bad_lexical_cast& e) {
        for (moduleLevel = static_cast<int>(ERROR);
             moduleLevel < NUM_LOG_LEVELS;
             ++moduleLevel)
        {
            if (level == logLevelNames[moduleLevel])
                break;
        }
        if (moduleLevel == NUM_LOG_LEVELS) {
            LOG(WARNING, "Ignoring bad log module level: %s", level.c_str());
            return;
        }
    }
    setLogLevel(static_cast<LogModule>(moduleIndex), moduleLevel);
#pragma GCC diagnostic pop
}

/**
 * Change the log level by a relative amount for a particular module.
 * \param[in] module
 *      The module whose level to change.
 * \param[in] delta
 *      The amount (positive or negative) by which to change the log level
 *      of \a module. The resulting level will be clamped to a valid LogLevel
 *      if it is out of range.
 */
void
Logger::changeLogLevel(LogModule module, int delta)
{
    Lock lock(mutex);
    LogLevel level = static_cast<LogLevel>(logLevels[module] + delta);
    logLevels[module] = level;
}

/**
 * Set the log level for all modules.
 * \param[in] level
 *      Messages for the modules at least as important as \a level will be
 *      logged.
 */
void
Logger::setLogLevels(LogLevel level)
{
    // No lock needed: doesn't access Logger object.
    for (int i = 0; i < NUM_LOG_MODULES; i++) {
        LogModule module = static_cast<LogModule>(i);
        setLogLevel(module, level);
    }
}
/**
 * Set the log level for all modules.
 * \param[in] level
 *      Messages for all modules at least as important as \a level will be
 *      logged. This will be clamped to a valid LogLevel if it is out of range.
 */
void
Logger::setLogLevels(int level)
{
    // No lock needed: doesn't access Logger object.
    if (level < 0)
        level = 0;
    else if (level >= NUM_LOG_LEVELS)
        level = NUM_LOG_LEVELS - 1;
    setLogLevels(static_cast<LogLevel>(level));
}
/**
 * Set the log level for all modules.
 * \param[in] level
 *      Messages for all modules at least as important as \a level will be
 *      logged. This will be clamped to a valid LogLevel if it is out of range.
 */
void
Logger::setLogLevels(string level)
{
// The uninitialized variable is inside boost::lexical_cast.
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
    // No lock needed: doesn't access Logger object.
    int moduleLevel;
    try {
        moduleLevel = boost::lexical_cast<int>(level);
    } catch (boost::bad_lexical_cast& e) {
        for (moduleLevel = static_cast<int>(ERROR);
             moduleLevel < NUM_LOG_LEVELS;
             ++moduleLevel)
        {
            if (level == logLevelNames[moduleLevel])
                break;
        }
        if (moduleLevel == NUM_LOG_LEVELS) {
            LOG(WARNING, "Ignoring bad log module level: %s", level.c_str());
            return;
        }
    }
    setLogLevels(moduleLevel);
#pragma GCC diagnostic pop
}

/**
 * Change the log level by a relative amount for all modules.
 * \param[in] delta
 *      The amount (positive or negative) by which to change the log level
 *      of the modules. The resulting level will be clamped to a valid LogLevel
 *      if it is out of range.
 */
void
Logger::changeLogLevels(int delta)
{
    // No lock needed: doesn't access Logger object.
    for (int i = 0; i < NUM_LOG_MODULES; i++) {
        LogModule module = static_cast<LogModule>(i);
        changeLogLevel(module, delta);
    }
}

/**
 * If there is currently a large amount of buffered log output that has
 * not yet been written, wait until this congestion eases before returning.
 * This method is intended for use by background loggers that generate
 * large amounts of data, such as the TimeTrace mechanism: they can
 * invoke this method occasionally to throttle themeselves so that there
 * is always plenty of buffer space for normal logging (if buffer space
 * runs out, normal logging will block, and this can lock up the server).
 */
void
Logger::waitIfCongested()
{
    Lock lock(mutex);
    while (1) {
        int bytesUsed = nextToInsert - nextToPrint;
        if (bytesUsed < 0) {
            bytesUsed = bytesUsed + bufferSize;
        }
        if (bytesUsed*2 < bufferSize) {
            return;
        }
        bufferSpaceFreed.wait(lock);
    }
}

/**
 * Log a backtrace for the system administrator.
 * This version doesn't provide C++ name demangling so it can be helpful
 * to run the output of the backtrace through c++filt.
 * \param[in] module
 *      The module to which the message pertains.
 * \param[in] level
 *      See #LOG.
 * \param[in] where
 *      The result of #HERE.
 */
void
Logger::logBacktrace(LogModule module, LogLevel level,
                     const CodeLocation& where)
{
    // No lock needed: doesn't access Logger object.
    const int maxFrames = 128;
    void* retAddrs[maxFrames];
    int frames = backtrace(retAddrs, maxFrames);
    char** symbols = backtrace_symbols(retAddrs, frames);
    if (symbols == NULL) {
        // If the malloc failed we might be able to get the backtrace out
        // to stderr still.
        backtrace_symbols_fd(retAddrs, frames, 2);
        return;
    }

    logMessage(false, module, level, where, "Backtrace:\n");
    for (int i = 0; i < frames; ++i)
        logMessage(false, module, level, where, "%s\n", symbols[i]);

    free(symbols);
}

/**
 * Log a message for the system administrator.
 * \param[in] collapse
 *      Collapse log messages when set.
 * \param[in] module
 *      The module to which the message pertains.
 * \param[in] level
 *      See #LOG.
 * \param[in] where
 *      The result of #HERE.
 * \param[in] fmt
 *      See #LOG except the string should end with a newline character.
 * \param[in] ...
 *      See #LOG.
 */
void
Logger::logMessage(bool collapse, LogModule module, LogLevel level,
                   const CodeLocation& where,
                   const char* fmt, ...)
{
    uint64_t start = Cycles::rdtsc();
    struct timespec now;
    clock_gettime(CLOCK_REALTIME, &now);
#ifdef TESTING
    if (testingLogTime != NULL) {
        now = *testingLogTime;
    }
#endif
    Lock lock(mutex);

    // See if this log message should be collapsed away entirely.
    SkipInfo* skip = NULL;
    int skipCount = 0;
    if (collapse) {
        auto messageId = std::make_pair(where.file, where.line);
        CollapseMap::iterator iter = collapseMap.find(messageId);
        if (iter != collapseMap.end()) {
            // We have printed this message before; if it was recent,
            // don't print the current message, but keep track of the fact
            // that we skipped it.
            skip = &iter->second;
            if (Util::timespecLess(now, skip->nextPrintTime)) {
                skip->skipCount++;
                return;
            }

            // We've printed this message before, but it was a while ago,
            // so print the message again.
            skipCount = skip->skipCount;
        } else {
            // Initialize a new entry in the map.
            if (collapseMap.size() >= maxCollapseMapSize) {
                // The map has gotten too large; just delete an entry
                // at random.
                collapseMap.erase(collapseMap.begin());
            }
            skip = &collapseMap[messageId];
        }
    }

#define MAX_MESSAGE_CHARS 2000
    // Extra space for a message about truncated characters, if needed.
#define TRUNC_MESSAGE_SPACE 50
    char buffer[MAX_MESSAGE_CHARS + TRUNC_MESSAGE_SPACE];
    int spaceLeft = MAX_MESSAGE_CHARS;
    if (testingBufferSize != 0) {
        spaceLeft = testingBufferSize;
    }
    int charsLost = 0;
    int charsWritten = 0;
    int actual;

    // Generate a message about discarded entries, if relevant.
    if (discardedEntries > 0) {
        CodeLocation here = HERE;
        actual = snprintf(buffer, spaceLeft,
                "%010lu.%09lu %s:%d in %s %s[%d]: %d log messages "
                "lost because of buffer overflow\n",
                now.tv_sec, now.tv_nsec, here.baseFileName(),
                here.line, here.function, logLevelNames[WARNING],
                ThreadId::get(), discardedEntries);
        if (actual >= spaceLeft) {
            // We ran out of space in the buffer (should never happen here).
            charsLost += 1 + actual - spaceLeft;
            actual = spaceLeft - 1;
        }
        charsWritten += actual;
        spaceLeft -= actual;
    }

    // Create the new log message in a local buffer. First write a standard
    // prefix containing timestamp, information about source file, etc.
    actual = snprintf(buffer+charsWritten, spaceLeft,
            "%010lu.%09lu %s:%d in %s %s[%d]: ",
            now.tv_sec, now.tv_nsec, where.baseFileName(), where.line,
            where.function, logLevelNames[level],
            ThreadId::get());
    if (actual >= spaceLeft) {
        // We ran out of space in the buffer (should never happen here).
        charsLost += 1 + actual - spaceLeft;
        actual = spaceLeft - 1;
    }
    charsWritten += actual;
    spaceLeft -= actual;

    // Next, add info about skipped entries, if relevant.
    if (skipCount > 0) {
        actual = snprintf(buffer + charsWritten, spaceLeft,
                "(%d duplicates of this message were skipped) ", skipCount);
        if (actual >= spaceLeft) {
            // We ran out of space in the buffer (should never happen here).
            charsLost += 1 + actual - spaceLeft;
            actual = spaceLeft - 1;
        }
        charsWritten += actual;
        spaceLeft -= actual;
    }

    // Last, add the caller's log message.
    va_list ap;
    va_start(ap, fmt);
    actual = vsnprintf(buffer + charsWritten, spaceLeft, fmt, ap);
    va_end(ap);
    if (actual >= spaceLeft) {
        // We ran out of space in the buffer.
        charsLost += 1 + actual - spaceLeft;
        actual = spaceLeft - 1;
    }
    charsWritten += actual;
    spaceLeft -= actual;

    if (charsLost > 0) {
        // Ran out of space: add a note about the lost info.
        charsWritten += snprintf(buffer + charsWritten, TRUNC_MESSAGE_SPACE,
                "... (%d chars truncated)\n", charsLost);
    }

    if (addToBuffer(buffer, charsWritten)) {
        discardedEntries = 0;
    } else {
        discardedEntries++;
    }

    if (collapse) {
        // If we reach this point of control, then we have just printed a
        // collapsible message (either it's new or  we haven't printed it in
        // a while). Save information so we don't print this message again
        // for a while.
        if (skip->message.empty()) {
            skip->message.assign(buffer, charsWritten);
        }
        skip->skipCount = 0;
        skip->nextPrintTime = Util::timespecAdd(now,
                {collapseIntervalMs/1000,
                (collapseIntervalMs%1000)*1000000});
        if (Util::timespecLess(skip->nextPrintTime, nextCleanTime)) {
            nextCleanTime = skip->nextPrintTime;
        }
    }

    // Print previously-deferred messages, if needed.
    if (Util::timespecLessEqual(nextCleanTime, now)) {
        cleanCollapseMap(now);
    }

    // Make sure this method did not take very long to execute.  If the
    // logging gets backed up it is really bad news, because it can lock
    // up the server so that it appears dead and crash recovery happens
    // (e.g., the Dispatch thread might be blocked waiting to log a message).
    // Thus, generate a log message to help people trying to debug the "crash".
    double elapsedMs = Cycles::toSeconds(Cycles::rdtsc() - start)*1e3;
    if (elapsedMs > 10) {
        struct timespec now;
        clock_gettime(CLOCK_REALTIME, &now);
        CodeLocation here = HERE;
        snprintf(buffer, sizeof(buffer), "%010lu.%09lu %s:%d in %s "
                "ERROR[%d]: Logger got stuck for %.1f ms, which could "
                "hang server\n",
                now.tv_sec, now.tv_nsec, here.baseFileName(),
                here.line, here.function, ThreadId::get(),
                elapsedMs);
        addToBuffer(buffer, downCast<int>(strlen(buffer)));
    }
}

/**
 * This method is called to print delayed log messages. If a bunch of
 * duplicates for a message were suppressed, but no new copies of that
 * message are logged, we eventually want to print out information about
 * the suppressed duplicates. This method does that (and it also removes
 * old entries from \c collapseMap). The caller must hold the Logger lock.
 *
 * \param now
 *      The current time: information will be printed about any suppressed
 *      duplicates with \c nextPrintTime less than or equal to this, or
 *      such entries will be removed from collapseMap if there are no
 *      suppressed duplicates.
 */
void
Logger::cleanCollapseMap(struct timespec now)
{
    // No lock needed: already locked by caller

    // We won't check again for a very long time, unless there are
    // entries in the collapse map (see below).
    nextCleanTime = Util::timespecAdd(now, {1000, 0});
    CollapseMap::iterator iter = collapseMap.begin();
    while (iter != collapseMap.end()) {
        auto messageId = &iter->first;
        SkipInfo* skip = &iter->second;
        if (Util::timespecLessEqual(skip->nextPrintTime, now)) {
            // This entry has been around for a while: if there were
            // suppressed messages then print one; otherwise delete the
            // entry.
            if (skip->skipCount == 0) {
                // This entry is old and there haven't been any suppressed
                // log messages for it; just delete the entry.
                iter++;
                collapseMap.erase(*messageId);
                continue;
            }

            // This entry contains suppressed log messages, but it's been
            // a long time since we printed the last one; print another one.
            // (splice a message about skipped messages into the log message).
            size_t i = skip->message.find("]:");
            assert(i > 20);
            string newMessage = format("%010lu.%09lu", now.tv_sec, now.tv_nsec)
                    + skip->message.substr(20, i+2-20)
                    + format(" (%d duplicates of this message were skipped)",
                    skip->skipCount)
                    + skip->message.substr(i+2);
            if (addToBuffer(newMessage.c_str(),
                    downCast<int>(newMessage.size()))) {
                skip->skipCount = 0;
            }
            skip->nextPrintTime = Util::timespecAdd(now,
                    {collapseIntervalMs/1000,
                    (collapseIntervalMs%1000)*1000000});
        }
        if (Util::timespecLess(skip->nextPrintTime, nextCleanTime)) {
            nextCleanTime = skip->nextPrintTime;
        }
        iter++;
    }
}

/**
 * This method copies a block of data into the printBuffer and wakes up
 * the print thread to output it.
 *
 * \param src
 *      First byte in block of data to add to the buffer.
 * \param length
 *      Total number of bytes of data to add.  Must be > 0.
 *
 * \return
 *      The return value is true if the data was successfully added, and false
 *      if there wasn't enough space for all of it (in which case none is
 *      added).
 */
bool
Logger::addToBuffer(const char* src, int length)
{
    // No lock needed: already locked by caller

    // First, write data at the end of the buffer, if there's space there.
    if (nextToInsert >= nextToPrint) {
        int count = std::min(bufferSize - nextToInsert, length);
        memcpy(messageBuffer + nextToInsert, src, count);
        src += count;
        nextToInsert += count;
        if (nextToInsert == bufferSize) {
            nextToInsert = 0;
        }
        length -= count;
    }

    if (length > 0) {
        // We get here if the space at the end of the buffer was all
        // occupied, or if that wasn't enough to hold all of the data.
        // See if there's space at the beginning of the buffer. Note: we
        // can't use the last byte before nextToPrint; if we do, we'll
        // wrap around and lose all the buffered data.
        if (length > (nextToPrint - 1 - nextToInsert)) {
            return false;
        }
        memcpy(messageBuffer + nextToInsert, src, length);
        nextToInsert += length;
    }

    if (!testingNoNotify) {
        logDataAvailable.notify_one();
    }
    return true;
}

/**
 * This method is the main program for a separate thread that runs in the
 * background to print log messages to secondary storage (this way, the
 * main RAMCloud threads aren't delayed for I/O).
 *
 * \param logger
 *      The owning Logger. This thread accesses only information related
 *      to the buffer and the output file.
 */
void
Logger::printThreadMain(Logger* logger)
{
    Lock lock(logger->mutex);
    int bytesToPrint;
    ssize_t bytesPrinted;
    while (true) {
        if (logger-> printThreadExit) {
            TEST_LOG("print thread exiting");
            return;
        }

        // Handle buffer wraparound.
        if (logger->nextToPrint >= logger->bufferSize) {
            logger->nextToPrint = 0;
        }

        // See if there is new log data to print.
        bytesToPrint = logger->nextToInsert - logger->nextToPrint;
        if (bytesToPrint < 0) {
            // If we get here, it means that nextToInsert has wrapped back
            // to the start of the buffer.
            bytesToPrint = logger->bufferSize - logger->nextToPrint;
        }
        if (bytesToPrint == 0) {
            // The line below is not needed for correctness, but it
            // results in better cache performance: as long as the printer
            // keeps up, only the first part of the buffer will be used,
            // and the later parts will never be touched.
            logger->nextToPrint = logger->nextToInsert = 0;
            logger->logDataAvailable.wait(lock);
            continue;
        }

        // Print whatever data is available (be careful not to hold the
        // buffer lock during expensive kernel calls, since this will
        // block calls to logMessage).
        {
            int fd = logger->fd;
            int nextToPrint = logger->nextToPrint;
            Unlock<SpinLock> unlock(logger->mutex);
            // Since we release the lock before write, there could be a race
            // condition when another thread acquires the lock and tries to
            // close fd. Therefore, we must ensure that setLogFile() will only
            // be called when the print thread is not running.
            bytesPrinted = write(fd, logger->messageBuffer + nextToPrint,
                    bytesToPrint);
            while (logger->testingStallPrintThread) {
                // Empty loop body.
            }
        }
        if (bytesPrinted < 0) {
            fprintf(stderr, "Error writing log: %s\n", strerror(errno));

            // Skip these bytes; otherwise we'll loop infinitely.
            logger->nextToPrint += downCast<int>(bytesToPrint);
        } else {
            logger->nextToPrint += downCast<int>(bytesPrinted);
        }
        logger->bufferSpaceFreed.notify_all();
    }
}

/**
 * Restore a logger to its default initialized state. Used primarily by tests.
 */
void
Logger::reset()
{
    Lock lock(mutex);
    if (mustCloseFd) {
        close(fd);
    }
    fd = 2;
    mustCloseFd = 0;
    collapseMap.clear();
    collapseIntervalMs = DEFAULT_COLLAPSE_INTERVAL;
    maxCollapseMapSize = DEFAULT_COLLAPSE_MAP_LIMIT;
    nextCleanTime = {0, 0};
    nextToInsert = 0;
    nextToPrint = 0;
    discardedEntries = 0;
    testingBufferSize = 0;
    testingLogTime = NULL;
    testingNoNotify = false;
}

/**
 * Wait for all buffered log messages to be printed. This method is intended
 * only for tests and a few special situations such as application exit. It
 * should *not* be used in the normal course of logging messages, since it
 * can result in long delays that could potentially cause the server to be
 * considered crashed.
 */
void
Logger::sync()
{
    Lock lock(mutex);
    while (nextToInsert != nextToPrint) {
        Unlock<SpinLock> unlock(mutex);
        usleep(100);
    }
}

/**
 * This replaces the default __assert_fail function, which is invoked by the
 * assert macro when there are errors. This function differs from the default
 * __assert_fail function in that it dumps its output to the RAMCloud log file.
 *
 * \param assertion
 *      Textual description of the assertion that failed.
 * \param file
 *      Name of the source file containing the assertion.
 * \param line
 *      Line number of the line containing the assertion, within \a file.
 * \param function
 *      Textual description of the function containing the assertion.
 */
void
Logger::assertionError(const char *assertion, const char *file,
                       unsigned int line, const char *function)
{
    DIE("Assertion `%s' failed at %s:%u in %s",
            assertion, file, line, function);
}

/**
 * Signal handler which logs a backtrace and exits.
 *
 * \param signal
 *      Signal number which caused the handler to be invoked.
 * \param info
 *      Details about the cause of the signal. Used to find the
 *      faulting address for segfaults.
 * \param ucontext
 *      CPU context at the time the signal occurred.
 */
static void
criticalErrorHandler(int signal, siginfo_t* info, void* ucontext)
{
    ucontext_t* uc = static_cast<ucontext_t*>(ucontext);
    void* callerAddress =
        reinterpret_cast<void*>(uc->uc_mcontext.gregs[REG_RIP]);

    LOG(ERROR, "Signal %d (%s) at address %p from %p",
        signal, strsignal(signal), info->si_addr,
        callerAddress);

    const int maxFrames = 128;
    void* retAddrs[maxFrames];
    int frames = backtrace(retAddrs, maxFrames);

    // Overwrite sigaction with caller's address.
    retAddrs[1] = callerAddress;

    char** symbols = backtrace_symbols(retAddrs, frames);
    if (symbols == NULL) {
        // If the malloc failed we might be able to get the backtrace out
        // to stderr still.
        backtrace_symbols_fd(retAddrs, frames, 2);
        return;
    }

    LOG(ERROR, "Backtrace:");
    for (int i = 1; i < frames; ++i)
        LOG(ERROR, "%s\n", symbols[i]);
    Logger::get().sync();

    free(symbols);

    // use abort, rather than exit, to dump core/trap in gdb
    abort();
}

/**
 * Logs a backtrace and exits. For use with std::set_terminate().
 */
static void
terminateHandler()
{
    BACKTRACE(ERROR);

    // use abort, rather than exit, to dump core/trap in gdb
    abort();
}

/**
 * Install handlers for SIGSEGV and std::terminate which log backtraces to
 * the process log file before exiting. This function currently ignores all
 * the thread-safety issues that signals introduce. Since the handlers are
 * only invoked in fatal cases its probably not a big deal.
 */
void
Logger::installCrashBacktraceHandlers()
{
    struct sigaction signalAction;
    signalAction.sa_sigaction = criticalErrorHandler;
    signalAction.sa_flags = SA_RESTART | SA_SIGINFO;

    if (sigaction(SIGSEGV, &signalAction, NULL) != 0)
        LOG(ERROR, "Couldn't set signal handler for SIGSEGV, oh well");

    std::set_terminate(terminateHandler);
}

} // end RAMCloud

/**
 * This replaces the default __assert_fail function, which is invoked by the
 * assert macro when there are errors. This function differs from the default
 * __assert_fail function in that it dumps its output to the RAMCloud log file.
 *
 * \param assertion
 *      Textual description of the assertion that failed.
 * \param file
 *      Name of the source file containing the assertion.
 * \param line
 *      Line number of the line containing the assertion, within \a file.
 * \param function
 *      Textual description of the function containing the assertion.
 */
void
__assert_fail(const char *assertion, const char *file, unsigned int line,
        const char *function)
{
    RAMCloud::Logger::get().assertionError(assertion, file, line, function);
    abort();             // To keep the compiler happy; never gets executed.
}

ext: cc

max_issues_repo_name: rodgarrison/RAMCloud

max_forks_repo_name: rodgarrison/RAMCloud


text:
#include <boost/smart_ptr/intrusive_ref_counter.hpp>
#include "IntrusivePtr.hpp"

struct Item : public boost::intrusive_ref_counter< Item, boost::thread_safe_counter > {
	int value;
	Item(int v = 0) : value(v) {}
};

typedef IntrusivePtr<Item> ItemPtr;

int main() {
	IntrusivePtr<Item> p1(new Item(12));
	IntrusivePtr<Item> p2(new Item(23));
	IntrusivePtr<Item> p3(new Item(34));

	p2 = NULL;
	p2 = p1;

	for(int i = 0; i < 100; ++i) {
		ItemPtr _tmp(p2);
	}

}

max_stars_repo_name: amikey/music-player-core

max_issues_repo_name: n-peugnet/music-player-core

max_forks_repo_name: n-peugnet/music-player-core


text:
/********************************************************************************
 * Copyright 2009 The Robotics Group, The Maersk Mc-Kinney Moller Institute,
 * Faculty of Engineering, University of Southern Denmark
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 ********************************************************************************/

#ifndef RW_MATH_INERTIAMATRIX_HPP
#define RW_MATH_INERTIAMATRIX_HPP

/**
 * @file InertiaMatrix.hpp
 */

#if !defined(SWIG)
#include "Rotation3D.hpp"
#include "Vector3D.hpp"

#include <rw/common/Serializable.hpp>

#include <Eigen/Core>
#endif

namespace rw { namespace math {

    /** @addtogroup math */
    /* @{*/

    /**
     * @brief A 3x3 inertia matrix
     */
    template< class T = double > class InertiaMatrix
    {
      public:
        //! @brief The type of the internal Eigen matrix implementation.
        typedef Eigen::Matrix< T, 3, 3 > Base;
#if !defined(SWIGJAVA)

        /**
         * @brief Constructs an initialized 3x3 rotation matrix
         *
         * @param r11 \f$ r_{11} \f$
         * @param r12 \f$ r_{12} \f$
         * @param r13 \f$ r_{13} \f$
         * @param r21 \f$ r_{21} \f$
         * @param r22 \f$ r_{22} \f$
         * @param r23 \f$ r_{23} \f$
         * @param r31 \f$ r_{31} \f$
         * @param r32 \f$ r_{32} \f$
         * @param r33 \f$ r_{33} \f$
         *
         * @f$
         *  \mathbf{R} =
         *  \left[
         *  \begin{array}{ccc}
         *  r_{11} & r_{12} & r_{13} \\
         *  r_{21} & r_{22} & r_{23} \\
         *  r_{31} & r_{32} & r_{33}
         *  \end{array}
         *  \right]
         * @f$
         */

#endif
        InertiaMatrix (T r11, T r12, T r13, T r21, T r22, T r23, T r31, T r32, T r33)
        {
            _matrix (0, 0) = r11;
            _matrix (0, 1) = r12;
            _matrix (0, 2) = r13;
            _matrix (1, 0) = r21;
            _matrix (1, 1) = r22;
            _matrix (1, 2) = r23;
            _matrix (2, 0) = r31;
            _matrix (2, 1) = r32;
            _matrix (2, 2) = r33;
        }
#if !defined(SWIGJAVA)

        /**
         * @brief Constructs an initialized 3x3 rotation matrix
         * @f$ \robabx{a}{b}{\mathbf{R}} =
         * \left[
         *  \begin{array}{ccc}
         *   \robabx{a}{b}{\mathbf{i}} & \robabx{a}{b}{\mathbf{j}} & \robabx{a}{b}{\mathbf{k}}
         *  \end{array}
         * \right]
         * @f$
         *
         * @param i @f$ \robabx{a}{b}{\mathbf{i}} @f$
         * @param j @f$ \robabx{a}{b}{\mathbf{j}} @f$
         * @param k @f$ \robabx{a}{b}{\mathbf{k}} @f$
         */

#endif
        InertiaMatrix (const rw::math::Vector3D< T >& i, const rw::math::Vector3D< T >& j,
                       const rw::math::Vector3D< T >& k)
        {
            _matrix (0, 0) = i[0];
            _matrix (0, 1) = j[0];
            _matrix (0, 2) = k[0];
            _matrix (1, 0) = i[1];
            _matrix (1, 1) = j[1];
            _matrix (1, 2) = k[1];
            _matrix (2, 0) = i[2];
            _matrix (2, 1) = j[2];
            _matrix (2, 2) = k[2];
        }

        /**
         * @brief constructor - where only the diagonal is set
         * @param i [in] m(0,0)
         * @param j [in] m(1,1)
         * @param k [in] m(2,2)
         */
        InertiaMatrix (T i = 0.0, T j = 0.0, T k = 0.0)
        {
            _matrix (0, 0) = i;
            _matrix (0, 1) = 0;
            _matrix (0, 2) = 0;
            _matrix (1, 0) = 0;
            _matrix (1, 1) = j;
            _matrix (1, 2) = 0;
            _matrix (2, 0) = 0;
            _matrix (2, 1) = 0;
            _matrix (2, 2) = k;
        }

        /**
         * @brief Construct an internal matrix from a Eigen::MatrixBase
         * It is the responsibility of the user that 3x3 matrix is indeed an
         * inertia matrix.
         */
        explicit InertiaMatrix (const Base& r) : _matrix (r) {}

        /**
         * @brief Returns reference to matrix element
         * @param row [in] row
         * @param column [in] column
         * @return reference to the element
         */
        T& operator() (size_t row, size_t column) { return _matrix (row, column); }

        /**
         * @brief Returns reference to matrix element
         * @param row [in] row
         * @param column [in] column
         * @return reference to the element
         */
        const T& operator() (size_t row, size_t column) const { return _matrix (row, column); }
#if defined(SWIG)
        MATRIXOPERATOR (T);
#endif
        /**
         * @brief Returns reference to the internal 3x3 matrix
         */
        const Base& e () const { return _matrix; }

        /**
         * @brief Returns reference to the internal 3x3 matrix
         */
        Base& e () { return _matrix; }

#if !defined(SWIG)
        /**
         * @brief Calculates \f$ \robabx{a}{c}{\mathbf{R}} =
         * \robabx{a}{b}{\mathbf{R}} \robabx{b}{c}{\mathbf{R}} \f$
         * @param aRb [in] \f$ \robabx{a}{b}{\mathbf{R}} \f$
         * @param bRc [in] \f$ \robabx{b}{c}{\mathbf{R}} \f$
         * @return \f$ \robabx{a}{c}{\mathbf{R}} \f$
         */
        friend InertiaMatrix operator* (const rw::math::Rotation3D< T >& aRb,
                                        const InertiaMatrix& bRc)
        {
            return InertiaMatrix (aRb.e () * bRc.e ());
        }
#endif
#if !defined(SWIGJAVA)

        /**
         * @brief Calculates \f$ \robabx{a}{c}{\mathbf{R}} =
         * \robabx{a}{b}{\mathbf{R}} \robabx{b}{c}{\mathbf{R}} \f$
         * @param bRc [in] \f$ \robabx{b}{c}{\mathbf{R}} \f$
         * @return \f$ \robabx{a}{c}{\mathbf{R}} \f$
         */
#endif
        InertiaMatrix operator* (const rw::math::Rotation3D< T >& bRc) const
        {
            return InertiaMatrix (this->e () * bRc.e ());
        }

        /**
         * @brief Calculates the addition between the two InertiaMatrices
         */
        InertiaMatrix operator+ (const InertiaMatrix& I2) const
        {
            return InertiaMatrix (this->e () + I2.e ());
        }

        /**
         * @brief Calculates \f$ \robabx{a}{c}{\mathbf{v}} =
         * \robabx{a}{b}{\mathbf{R}} \robabx{b}{c}{\mathbf{v}} \f$
         * @param bVc [in] \f$ \robabx{b}{c}{\mathbf{v}} \f$
         * @return \f$ \robabx{a}{c}{\mathbf{v}} \f$
         */
        rw::math::Vector3D< T > operator* (const rw::math::Vector3D< T >& bVc) const
        {
            return rw::math::Vector3D< T > (this->e () * bVc.e ());
        }
#if !defined(SWIG)
        /**
         * @brief Writes rotation matrix to stream
         * @param os [in/out] output stream to use
         * @param r [in] rotation matrix to print
         * @return the updated output stream
         */
        friend std::ostream& operator<< (std::ostream& os, const InertiaMatrix& r)
        {
            return os << r.e ();
        }
#else
        TOSTRING (rw::math::InertiaMatrix< T >);
#endif

        /**
         * @brief Make inertia matrix for a solid sphere.
         * @param mass [in] mass of solid sphere.
         * @param radi [in] radius of sphere.
         * @return the inertia matrix.
         */
        static InertiaMatrix< T > makeSolidSphereInertia (T mass, T radi)
        {
            T tmpV = (T) (2.0 / 5.0) * mass * radi * radi;
            return InertiaMatrix< T > (tmpV, 0, 0, 0, tmpV, 0, 0, 0, tmpV);
        }

        /**
         * @brief Make inertia matrix for a hollow sphere.
         * @param mass [in] mass of hollow sphere.
         * @param radi [in] radius of sphere.
         * @return the inertia matrix.
         */
        static InertiaMatrix< T > makeHollowSphereInertia (T mass, T radi)
        {
            T tmpV = (T) (2.0 / 3.0) * mass * radi * radi;
            return InertiaMatrix< T > (tmpV, 0, 0, 0, tmpV, 0, 0, 0, tmpV);
        }

        /**
         * @brief calculates the inertia of a cuboid where the reference frame is in the
         * center of the cuboid with
         * @param mass
         * @param x
         * @param y
         * @param z
         * @return
         */
        static InertiaMatrix< T > makeCuboidInertia (T mass, T x, T y, T z)
        {
            return InertiaMatrix< T > ((T) (1 / 12.0 * mass * (y * y + z * z)),
                                       0,
                                       0,
                                       0,
                                       (T) (1 / 12.0 * mass * (x * x + z * z)),
                                       0,
                                       0,
                                       0,
                                       (T) (1 / 12.0 * mass * (x * x + y * y)));
        }

      private:
        Base _matrix;
    };
#if !defined(SWIGJAVA)

    /**
     * @brief Calculates the inverse @f$ \robabx{b}{a}{\mathbf{R}} =
     * \robabx{a}{b}{\mathbf{R}}^{-1} @f$ of a rotation matrix
     *
     * @param aRb [in] the rotation matrix @f$ \robabx{a}{b}{\mathbf{R}} @f$
     *
     * @return the matrix inverse @f$ \robabx{b}{a}{\mathbf{R}} =
     * \robabx{a}{b}{\mathbf{R}}^{-1} @f$
     *
     * @f$ \robabx{b}{a}{\mathbf{R}} = \robabx{a}{b}{\mathbf{R}}^{-1} =
     * \robabx{a}{b}{\mathbf{R}}^T @f$
     */

    #endif
    template< class Q > InertiaMatrix< Q > inverse (const InertiaMatrix< Q >& aRb)
    {
        return InertiaMatrix< Q > (aRb.e ().inverse ());
    }

    /**
     * @brief Casts InertiaMatrix<T> to InertiaMatrix<Q>
     * @param rot [in] InertiaMatrix with type T
     * @return InertiaMatrix with type Q
     */
    template< class Q, class T > InertiaMatrix< Q > cast (const InertiaMatrix< T >& rot)
    {
        InertiaMatrix< Q > res;
        for (size_t i = 0; i < 3; i++)
            for (size_t j = 0; j < 3; j++)
                res (i, j) = static_cast< Q > (rot (i, j));
        return res;
    }
#if !defined(SWIG)
    extern template class rw::math::InertiaMatrix< double >;
    extern template class rw::math::InertiaMatrix< float >;
#else
    SWIG_DECLARE_TEMPLATE (InertiaMatrixd, rw::math::InertiaMatrix< double >);
    SWIG_DECLARE_TEMPLATE (InertiaMatrixf, rw::math::InertiaMatrix< float >);
#endif

    using InertiaMatrixd = InertiaMatrix< double >;
    using InertiaMatrixf = InertiaMatrix< float >;

    /*@}*/
}}    // namespace rw::math

namespace rw { namespace common {
    class OutputArchive;
    class InputArchive;
    namespace serialization {
        /**
         * @copydoc rw::common::serialization::write
         * @relatedalso rw::math::InertiaMatrix
         */
        template<>
        void write (const rw::math::InertiaMatrix< double >& sobject,
                    rw::common::OutputArchive& oarchive, const std::string& id);

        /**
         * @copydoc rw::common::serialization::write
         * @relatedalso rw::math::InertiaMatrix
         */
        template<>
        void write (const rw::math::InertiaMatrix< float >& sobject,
                    rw::common::OutputArchive& oarchive, const std::string& id);

        /**
         * @copydoc rw::common::serialization::read
         * @relatedalso rw::math::InertiaMatrix
         */
        template<>
        void read (rw::math::InertiaMatrix< double >& sobject, rw::common::InputArchive& iarchive,
                   const std::string& id);

        /**
         * @copydoc rw::common::serialization::read
         * @relatedalso rw::math::InertiaMatrix
         */
        template<>
        void read (rw::math::InertiaMatrix< float >& sobject, rw::common::InputArchive& iarchive,
                   const std::string& id);
    }    // namespace serialization
}}       // namespace rw::common

#endif    // end include guard


text:
/***************************************************************************
 * Copyright 1998-2018 by authors (see AUTHORS.txt)                        *
 *                                                                         *
 *   This file is part of LuxCoreRender.                                   *
 *                                                                         *
 * Licensed under the Apache License, Version 2.0 (the "License");         *
 * you may not use this file except in compliance with the License.        *
 * You may obtain a copy of the License at                                 *
 *                                                                         *
 *     http://www.apache.org/licenses/LICENSE-2.0                          *
 *                                                                         *
 * Unless required by applicable law or agreed to in writing, software     *
 * distributed under the License is distributed on an "AS IS" BASIS,       *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.*
 * See the License for the specific language governing permissions and     *
 * limitations under the License.                                          *
 ***************************************************************************/

#include <boost/lexical_cast.hpp>

#include "luxrays/core/color/color.h"
#include "slg/samplers/sampler.h"
#include "slg/samplers/metropolis.h"

using namespace std;
using namespace luxrays;
using namespace slg;

//------------------------------------------------------------------------------
// MetropolisSamplerSharedData
//------------------------------------------------------------------------------

MetropolisSamplerSharedData::MetropolisSamplerSharedData() : SamplerSharedData() {
	totalLuminance = 0.;
	sampleCount = 0.;
}

SamplerSharedData *MetropolisSamplerSharedData::FromProperties(const Properties &cfg,
		RandomGenerator *rndGen, Film *film) {
	return new MetropolisSamplerSharedData();
}

//------------------------------------------------------------------------------
// Metropolis sampler
//------------------------------------------------------------------------------

MetropolisSampler::MetropolisSampler(RandomGenerator *rnd, Film *flm,
		const FilmSampleSplatter *flmSplatter, const u_int maxRej,
		const float pLarge, const float imgRange,
		MetropolisSamplerSharedData *samplerSharedData) : Sampler(rnd, flm, flmSplatter),
		sharedData(samplerSharedData),
		maxRejects(maxRej),	largeMutationProbability(pLarge), imageMutationRange(imgRange),
		samples(NULL), sampleStamps(NULL), currentSamples(NULL), currentSampleStamps(NULL),
		cooldown(true) {
}

MetropolisSampler::~MetropolisSampler() {
	delete[] samples;
	delete[] sampleStamps;
	delete[] currentSamples;
	delete[] currentSampleStamps;
}

// Mutate a value in the range [0-1]
static float Mutate(const float x, const float randomValue) {
	static const float s1 = 1.f / 512.f;
	static const float s2 = 1.f / 16.f;

	const float dx = s1 / (s1 / s2 + fabsf(2.f * randomValue - 1.f)) -
			s1 / (s1 / s2 + 1.f);

	float mutatedX = x;
	if (randomValue < .5f) {
		mutatedX += dx;
		mutatedX = (mutatedX < 1.f) ? mutatedX : (mutatedX - 1.f);
	} else {
		mutatedX -= dx;
		mutatedX = (mutatedX < 0.f) ? (mutatedX + 1.f) : mutatedX;
	}

	// mutatedX can still be 1.f due to numerical precision problems
	if (mutatedX == 1.f)
		mutatedX = 0.f;

	return mutatedX;
}

// Mutate a value max. by a range value
float MutateScaled(const float x, const float range, const float randomValue) {
	static const float s1 = 32.f;

	const float dx = range / (s1 / (1.f + s1) + (s1 * s1) / (1.f + s1) *
		fabs(2.f * randomValue - 1.f)) - range / s1;

	float mutatedX = x;
	if (randomValue < .5f) {
		mutatedX += dx;
		mutatedX = (mutatedX < 1.f) ? mutatedX : (mutatedX - 1.f);
	} else {
		mutatedX -= dx;
		mutatedX = (mutatedX < 0.f) ? (mutatedX + 1.f) : mutatedX;
	}

	// mutatedX can still be 1.f due to numerical precision problems
	if (mutatedX == 1.f)
		mutatedX = 0.f;

	return mutatedX;
}

void MetropolisSampler::RequestSamples(const u_int size) {
	sampleSize = size;
	samples = new float[sampleSize];
	sampleStamps = new u_int[sampleSize];
	currentSamples = new float[sampleSize];
	currentSampleStamps = new u_int[sampleSize];

	isLargeMutation = true;
	weight = 0.f;
	consecRejects = 0;
	currentLuminance = 0.f;
	fill(sampleStamps, sampleStamps + sampleSize, 0);
	stamp = 1;
	currentStamp = 1;
	currentSampleResult.resize(0);
}

float MetropolisSampler::GetSample(const u_int index) {
	assert (index < sampleSize);

	u_int sampleStamp = sampleStamps[index];

	float s;
	if (sampleStamp == 0) {
		s = rndGen->floatValue();
		sampleStamp = 1;

		assert (s != 1.f);
	} else
		s = samples[index];

	// Mutate the sample up to the currentStamp
	if (film && ((index == 0) || (index == 1))) {
		// 0 and 1 are used for image X/Y
		for (u_int i = sampleStamp; i < stamp; ++i)
			s = MutateScaled(s, imageMutationRange, rndGen->floatValue());

		samples[index] = s;
		sampleStamps[index] = stamp;

		const u_int *subRegion = film->GetSubRegion();
		const u_int subRegionIndex = (index == 0) ? 0 : 2;
		s = subRegion[subRegionIndex] + s * (subRegion[subRegionIndex + 1] - subRegion[subRegionIndex] + 1);

		return s;
	} else {
		for (u_int i = sampleStamp; i < stamp; ++i)
			s = Mutate(s, rndGen->floatValue());

		samples[index] = s;
		sampleStamps[index] = stamp;

		return s;
	}
}

void MetropolisSampler::NextSample(const vector<SampleResult> &sampleResults) {
	if (film)
		film->AddSampleCount(1.0);

	// Calculate the sample result luminance
	float newLuminance = 0.f;
	for (vector<SampleResult>::const_iterator sr = sampleResults.begin(); sr != sampleResults.end(); ++sr) {
		if (sr->HasChannel(Film::RADIANCE_PER_PIXEL_NORMALIZED)) {
			for (u_int i = 0; i < sr->radiance.size(); ++i) {
				const float luminance = sr->radiance[i].Y();
				assert (!isnan(luminance) && !isinf(luminance) && (luminance >= 0.f));

				if ((luminance > 0.f) && !isnan(luminance) && !isinf(luminance))
					newLuminance += luminance;
			}
		}

		if (sr->HasChannel(Film::RADIANCE_PER_SCREEN_NORMALIZED)) {
			for (u_int i = 0; i < sr->radiance.size(); ++i) {
				const float luminance = sr->radiance[i].Y();
				assert (!isnan(luminance) && !isinf(luminance) && (luminance >= 0.f));

				if ((luminance > 0.f) && !isnan(luminance) && !isinf(luminance))
					newLuminance += luminance;
			}
		}
	}

	if (isLargeMutation) {
		// Atomic 64bit for double are not available on all CPUs so I simply
		// avoid synchronization.
		sharedData->totalLuminance += newLuminance;
		sharedData->sampleCount += 1.;
	}

	const float invMeanIntensity = (sharedData->totalLuminance > 0.) ?
		static_cast<float>(sharedData->sampleCount / sharedData->totalLuminance) : 1.f;

	// Define the probability of large mutations. It is 50% if we are still
	// inside the cooldown phase.
	const float currentLargeMutationProbability = (cooldown) ? .5f : largeMutationProbability;

	// Calculate accept probability from old and new image sample
	float accProb;
	if ((currentLuminance > 0.f) && (consecRejects < maxRejects))
		accProb = Min<float>(1.f, newLuminance / currentLuminance);
	else
		accProb = 1.f;
	const float newWeight = accProb + (isLargeMutation ? 1.f : 0.f);
	weight += 1.f - accProb;

	const float rndVal = rndGen->floatValue();

	/*if (!cooldown && (currentSampleResult.size() > 0) && (sampleResults.size() > 0))
		printf("[%d] Current: (%f, %f, %f) [%f] Proposed: (%f, %f, %f) [%f] accProb: %f <%f>\n",
				consecRejects,
				currentSampleResult[0].radiance[0].c[0], currentSampleResult[0].radiance[0].c[1], currentSampleResult[0].radiance[0].c[2], weight,
				sampleResults[0].radiance[0].c[0], sampleResults[0].radiance[0].c[1], sampleResults[0].radiance[0].c[2], newWeight,
				accProb, rndVal);*/

	// Try or force accepting of the new sample
	if ((accProb == 1.f) || (rndVal < accProb)) {
		/*if (!cooldown)
			printf("\t\tACCEPTED !\n");*/

		// Add accumulated SampleResult of previous reference sample
		const float norm = weight / (currentLuminance * invMeanIntensity + currentLargeMutationProbability);
		if (norm > 0.f) {
			/*if (!cooldown)
				printf("\t\tContrib: (%f, %f, %f) [%f] consecutiveRejects: %d\n",
						currentSampleResult[0].radiance[0].c[0],
						currentSampleResult[0].radiance[0].c[1],
						currentSampleResult[0].radiance[0].c[2],
						norm, consecRejects);*/

			if (film)
				AtomicAddSamplesToFilm(currentSampleResult, norm);
		}

		lastSampleAcceptance = METRO_ACCEPTED;
		lastSampleWeight = norm;

		// Save new contributions for reference
		weight = newWeight;
		currentStamp = stamp;
		currentLuminance = newLuminance;
		copy(samples, samples + sampleSize, currentSamples);
		copy(sampleStamps, sampleStamps + sampleSize, currentSampleStamps);
		currentSampleResult = sampleResults;

		consecRejects = 0;
	} else {
		/*if (!cooldown)
			printf("\t\tREJECTED !\n");*/

		// Add contribution of new sample before rejecting it
		const float norm = newWeight / (newLuminance * invMeanIntensity + currentLargeMutationProbability);
		if (norm > 0.f) {
			/*if (!cooldown)
				printf("\t\tContrib: (%f, %f, %f) [%f] consecutiveRejects: %d\n",
						sampleResults[0].radiance[0].c[0],
						sampleResults[0].radiance[0].c[1],
						sampleResults[0].radiance[0].c[2],
						norm, consecRejects);*/

			if (film)
				AtomicAddSamplesToFilm(sampleResults, norm);
		}

		lastSampleAcceptance = METRO_REJECTED;
		lastSampleWeight = norm;

		// Restart from previous reference
		stamp = currentStamp;
		copy(currentSamples, currentSamples + sampleSize, samples);
		copy(currentSampleStamps, currentSampleStamps + sampleSize, sampleStamps);

		++consecRejects;
	}

	// Cooldown is used in order to not have problems in the estimation of meanIntensity
	// when large mutation probability is very small.
	if (cooldown) {
		// Check if it is time to end the cooldown (i.e. I have an average of
		// 1 sample for each pixel).
		const u_int pixelCount = film ? (film->GetWidth() * film->GetHeight()) : 8192;
		if (sharedData->sampleCount > pixelCount) {
			cooldown = false;
			isLargeMutation = (rndGen->floatValue() < currentLargeMutationProbability);
		} else
			isLargeMutation = (rndGen->floatValue() < .5f);
	} else
		isLargeMutation = (rndGen->floatValue() < currentLargeMutationProbability);

	if (isLargeMutation) {
		stamp = 1;
		fill(sampleStamps, sampleStamps + sampleSize, 0);
	} else
		++stamp;
}

// Used, most of the times, when not having a film
MetropolisSampleType MetropolisSampler::GetLastSampleAcceptance(float &weight) const {
	weight = lastSampleWeight;

	return lastSampleAcceptance;
}

Properties MetropolisSampler::ToProperties() const {
	return Sampler::ToProperties() <<
			Property("sampler.metropolis.largesteprate")(largeMutationProbability) <<
			Property("sampler.metropolis.maxconsecutivereject")(maxRejects) <<
			Property("sampler.metropolis.imagemutationrate")(imageMutationRange);
}

//------------------------------------------------------------------------------
// Static methods used by SamplerRegistry
//------------------------------------------------------------------------------

Properties MetropolisSampler::ToProperties(const Properties &cfg) {
	return Properties() <<
			cfg.Get(GetDefaultProps().Get("sampler.type")) <<
			cfg.Get(GetDefaultProps().Get("sampler.metropolis.largesteprate")) <<
			cfg.Get(GetDefaultProps().Get("sampler.metropolis.maxconsecutivereject")) <<
			cfg.Get(GetDefaultProps().Get("sampler.metropolis.imagemutationrate"));
}

Sampler *MetropolisSampler::FromProperties(const Properties &cfg, RandomGenerator *rndGen,
		Film *film, const FilmSampleSplatter *flmSplatter, SamplerSharedData *sharedData) {
	const float rate = Clamp(cfg.Get(GetDefaultProps().Get("sampler.metropolis.largesteprate")).Get<float>(), 0.f, 1.f);
	const u_int reject = cfg.Get(GetDefaultProps().Get("sampler.metropolis.maxconsecutivereject")).Get<u_int>();
	const float mutationRate = Clamp(cfg.Get(GetDefaultProps().Get("sampler.metropolis.imagemutationrate")).Get<float>(), 0.f, 1.f);

	return new MetropolisSampler(rndGen, film, flmSplatter,
			reject, rate, mutationRate,
			(MetropolisSamplerSharedData *)sharedData);
}

slg::ocl::Sampler *MetropolisSampler::FromPropertiesOCL(const Properties &cfg) {
	slg::ocl::Sampler *oclSampler = new slg::ocl::Sampler();

	oclSampler->type = slg::ocl::METROPOLIS;
	oclSampler->metropolis.largeMutationProbability = cfg.Get(GetDefaultProps().Get("sampler.metropolis.largesteprate")).Get<float>();
	oclSampler->metropolis.imageMutationRange = cfg.Get(GetDefaultProps().Get("sampler.metropolis.imagemutationrate")).Get<float>();
	oclSampler->metropolis.maxRejects = cfg.Get(GetDefaultProps().Get("sampler.metropolis.maxconsecutivereject")).Get<u_int>();

	return oclSampler;
}

Film::FilmChannelType MetropolisSampler::GetRequiredChannels(const luxrays::Properties &cfg) {
	return Film::NONE;
}

const Properties &MetropolisSampler::GetDefaultProps() {
	static Properties props = Properties() <<
			Sampler::GetDefaultProps() <<
			Property("sampler.type")(GetObjectTag()) <<
			Property("sampler.metropolis.largesteprate")(.4f) <<
			Property("sampler.metropolis.maxconsecutivereject")(512) <<
			Property("sampler.metropolis.imagemutationrate")(.1f);

	return props;
}

max_stars_repo_name: bartoszek/LuxCore

max_issues_repo_name: MetaRabbit/LuxCore

max_forks_repo_name: MetaRabbit/LuxCore


text:
#ifndef TREEDAG_DECOMPOSITIONDAG_HPP
#define TREEDAG_DECOMPOSITIONDAG_HPP

#include "separatorConfig.hpp"
#include "separation.hpp"

#include <boost/graph/adjacency_list.hpp>
#include <boost/unordered_map.hpp>
#include <boost/bimap.hpp>
#include <boost/bimap/unordered_set_of.hpp>

namespace treeDAG {

struct SeparatorNodeData : public SeparatorConfig
{
    static VertexIndexType NoInactiveComponentNumber() { return std::numeric_limits<VertexIndexType>::max(); }

    VertexSet separator;
    VertexSet inactiveComponents;
};

struct SubgraphNodeData : public SeparatorConfig
{
    VertexSet activeVertices;
    VertexSet otherVertices;
};

std::size_t hash_value(const SeparatorNodeData & separatorNode);
std::size_t hash_value(const SubgraphNodeData & subgraphNode);

bool operator==(const SeparatorNodeData & lhs, const SeparatorNodeData & rhs);
bool operator==(const SubgraphNodeData & lhs, const SubgraphNodeData & rhs);

std::ostream & operator<<(std::ostream & str, const SeparatorNodeData & separatorNode);
std::ostream & operator<<(std::ostream & str, const SubgraphNodeData & separatorNode);

struct DecompositionDAGNodeStreamWriter;

class DecompositionDAG : public boost::noncopyable, public SeparatorConfig
{
public:
    enum NodeType
    {
        NODE_Separator,
        NODE_Subgraph,
        NODE_Clique

    };

    typedef boost::adjacency_list<boost::hash_setS, boost::listS, boost::bidirectionalS, NodeType>                              Structure;
    typedef boost::graph_traits<Structure>::vertex_descriptor                                                               NodeDescriptor;

    // output methods
    void write_dot(std::ostream & stream) const;

    // structure methods
    const Structure & structure() const;
    NodeType nodeType(NodeDescriptor node) const;

    // find methods
    NodeDescriptor findSeparatorNode(const SeparatorNodeData & separatorNodeData) const;
    NodeDescriptor findSubgraphNode(const SubgraphNodeData & subgraphNodeData) const;
    const SeparatorNodeData * separatorNodeData(NodeDescriptor separatorNode) const;
    const SubgraphNodeData * subgraphNodeData(NodeDescriptor subgraphNode) const;

    // addition methods
    template <typename SubgraphNodeDataIterator>
    NodeDescriptor addSeparator(const SeparatorNodeData & separatorNodeData, SubgraphNodeDataIterator first, SubgraphNodeDataIterator last);
    NodeDescriptor addSubgraph(const SubgraphNodeData & subgraphNodeData);

    template <typename SeparatorNodeIterator>
    void addClique(NodeDescriptor subgraphNode, const VertexSet & clique, SeparatorNodeIterator first, SeparatorNodeIterator last);

    void cleanUp();

    std::size_t numberOfNodes() const { return boost::num_vertices(dag_); }
    std::size_t numberOfBranches() const { return boost::num_edges(dag_); }

    DecompositionDAGNodeStreamWriter nodeWriter(NodeDescriptor node) const;

    static NodeDescriptor InvalidNode() { return boost::graph_traits<Structure>::null_vertex(); }

private:
    friend class DecompositionDAGNodeStreamWriter;

    typedef boost::bimap<boost::bimaps::unordered_set_of<NodeDescriptor>, boost::bimaps::unordered_set_of<SubgraphNodeData > >  SubgraphMap;
    typedef boost::bimap<boost::bimaps::unordered_set_of<NodeDescriptor>, boost::bimaps::unordered_set_of<SeparatorNodeData > > SeparatorMap;
    typedef boost::unordered_map<NodeDescriptor, VertexSet>                                                                     CliqueSizeMap;

    void writeVertexSet(std::ostream & stream, const VertexSet & vertexSet) const;

    NodeDescriptor findOrCreateSeparatorNode(const SeparatorNodeData & separatorNode);
    NodeDescriptor findOrCreateSubgraphNode(const SubgraphNodeData & subgraphNode);

    void checkCliqueNode(NodeDescriptor node, boost::unordered_map<NodeDescriptor, std::size_t> & countMap);
    void checkSeparatorNode(NodeDescriptor node, boost::unordered_map<NodeDescriptor, std::size_t> & countMap);
    void checkSubgraphNode(NodeDescriptor node, boost::unordered_map<NodeDescriptor, std::size_t> & countMap);
    void cleanSubtree(NodeDescriptor node);
    void cleanupParallelEdges();

    bool hasSmallerTreewidthParallelEdge(NodeDescriptor subgraphNode, NodeDescriptor cliqueNode) const;

    Structure dag_;
    SubgraphMap subgraphMap_;
    SeparatorMap separatorMap_;
    CliqueSizeMap cliqueMap_;
};

struct DecompositionDAGNodeStreamWriter
{
    DecompositionDAGNodeStreamWriter() : dag_(0) {}

    friend std::ostream & operator<<(std::ostream & stream, const DecompositionDAGNodeStreamWriter & writer);

private:
    void toStream(std::ostream & stream) const;

    friend class DecompositionDAG;

    const DecompositionDAG * dag_;
    DecompositionDAG::NodeDescriptor node_;
};

} // namespace treeDAG

#include "decompositionDAG.hxx"

#endif // TREEDAG_DECOMPOSITIONDAG_HPP

max_stars_repo_name: thomasfannes/treeDAG

max_issues_repo_name: thomasfannes/treeDAG

max_forks_repo_name: thomasfannes/treeDAG