From 7970c0d5f2736ea973369c6bd6c4a919a058ce54 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sun, 17 Jan 2021 23:33:44 +0100 Subject: [PATCH 01/56] WIP - Finally working MLT implementation (much work still to do) --- src/spps-agh/calculation_cores/MLTCore.cpp | 104 ++++++++++++------ .../NextEventEstimationCore.cpp | 8 +- .../NextEventEstimationCore.h | 2 +- src/spps-agh/sppsNeeAGHTypes.h | 18 ++- src/spps-agh/tools/dotreflectionAGH.h | 24 ++-- 5 files changed, 104 insertions(+), 52 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index b1c6f7b6..382b662b 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -28,13 +28,13 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) CONF_PARTICULE_MLT inputParticle(configurationP); inputParticle.SetInitialState(inputParticle); - inputParticle.weight = 1/mutation_number; + inputParticle.weight = 0; int totalParticleNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); RunInitialSeed(inputParticle); Seeds.push_back(inputParticle); - b += inputParticle.totalProbability / totalParticleNumber; + b += inputParticle.probability.totalProbability() / totalParticleNumber; } else { @@ -65,6 +65,7 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) shadowRay = input_particle.shadowRays.front(); input_particle.shadowRays.pop_front(); //applicationTools.outputTool->NewParticule(confPart); + //shadowRay.energie *= input_particle.weight; shadowRay.energie *= input_particle.weight; Run(shadowRay); //applicationTools.outputTool->SaveParticule(); @@ -77,35 +78,43 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu double ratio; int mutation_counter = 0; + if (mutation_number == 0) + inputParticle.weight = b / inputParticle.probability.totalProbability(); + while (mutation_counter < mutation_number) { Mutate(inputParticle, mutatedParticle); Follow(mutatedParticle); - ratio = std::min(mutatedParticle.totalProbability / inputParticle.totalProbability, 1.); + ratio = std::min(mutatedParticle.probability.totalProbability() / inputParticle.probability.totalProbability(), 1.); - inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); - mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); + //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); + //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); + + // Rejection sampling + //inputParticle.weight += b / mutation_number / inputParticle.probability.totalProbability() * (1- ratio); + //mutatedParticle.weight += b / mutation_number / mutatedParticle.probability.totalProbability() * ratio; - //Cosnider if weight schould be different for each reflectio order! + // No rejection sampling + inputParticle.weight = b / inputParticle.probability.totalProbability(); + mutatedParticle.weight = b / mutatedParticle.probability.totalProbability(); if(mutatedParticle.largeStep==1) { - CastShadowRays(inputParticle); + //CastShadowRays(inputParticle); inputParticle = mutatedParticle; inputParticle.largeStep = 0; } else if(ratio > GetRandValue()) { - CastShadowRays(inputParticle); + //CastShadowRays(inputParticle); inputParticle = mutatedParticle; }else { - CastShadowRays(mutatedParticle); + //CastShadowRays(mutatedParticle); } mutation_counter++; } - CastShadowRays(inputParticle); } @@ -157,6 +166,8 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML { mutetedParticle = input_particle; mutetedParticle.shadowRays.clear(); + mutetedParticle.probability.partialProbability.clear(); + mutetedParticle.weight = 0; //bool large_step = (GetRandValue() < pLarge) ? 1 : 0; if (GetRandValue() < pLarge) @@ -194,14 +205,15 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn decimal absorption_atmospheric = configurationTool->freqList[configurationP.frequenceIndex]->absorption_atmospherique; decimal travelProbability = expf(-absorption_atmospheric * distance); - if (rndTravelProbability >= travelProbability) - { - configurationP.energie = 0; // La particule est détruite - configurationP.stateParticule = PARTICULE_STATE_ABS_ATMO; - configurationP.totalProbability *= (1 - travelProbability); - return true; - } - configurationP.totalProbability *= travelProbability; + //if (rndTravelProbability >= travelProbability) + //{ + // configurationP.energie = 0; // La particule est détruite + // configurationP.stateParticule = PARTICULE_STATE_ABS_ATMO; + // configurationP.probability.currentProbability *= (1 - travelProbability); + // return true; + //} + configurationP.energie *= travelProbability; + //configurationP.probability.currentProbability *= travelProbability; configurationP.totalDistance += distance; @@ -215,17 +227,18 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn t_Material_BFreq* materialInfo = &(*faceInfo).faceMaterial->matSpectrumProperty[configurationP.frequenceIndex]; //Test d'absorption en aléatoire - if (rndMatAbsorbtion <= materialInfo->absorption) - { - //Particule absorbée - if (configurationP.stateParticule == PARTICULE_STATE_ALIVE) - configurationP.stateParticule = PARTICULE_STATE_ABS_SURF; - configurationP.energie = 0.; - configurationP.totalProbability *= materialInfo->absorption; - return true; - } - - configurationP.totalProbability *= (1 - materialInfo->absorption); + //if (rndMatAbsorbtion <= materialInfo->absorption) + //{ + // //Particule absorbée + // if (configurationP.stateParticule == PARTICULE_STATE_ALIVE) + // configurationP.stateParticule = PARTICULE_STATE_ABS_SURF; + // configurationP.energie = 0.; + // configurationP.probability.currentProbability *= materialInfo->absorption; + // return true; + //} + + configurationP.energie *= (1 - materialInfo->absorption); + //configurationP.probability.currentProbability *= (1 - materialInfo->absorption); //******************************* UPDATE TIME STEP NUMBER ***********************************// double timestepNum = configurationP.totalDistance / configurationP.direction.length(); @@ -248,8 +261,15 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn else faceNormal = faceInfo->normal; + std::vector srProbabilities; + GenerateShadowRays(configurationP, materialInfo, faceNormal, deltaT, distanceToTravel, configurationP.shadowRays, faceInfo, &srProbabilities); - GenerateShadowRays(configurationP, materialInfo, faceNormal, deltaT, distanceToTravel, configurationP.shadowRays, faceInfo); + if (configurationP.reflectionOrder > 0) + { + for (double srProb : srProbabilities) { + configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability * srProb); + } + } //******************************* REFLECTION ***************************************// @@ -258,20 +278,33 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn double dirProbability; //Get direction for diffuse or specular part based on material info - if (materialInfo->diffusion == 1 || pSpecular*((1 - rndMatAbsorbtion) / materialInfo->absorption) < materialInfo->diffusion) + pSpecular = 1; + if (materialInfo->diffusion == 1 || pSpecular*((1 - rndMatAbsorbtion)) < materialInfo->diffusion) { nouvDirection = ReflectionLawsMLT::SolveDiffusePart(configurationP.direction, *materialInfo, faceNormal, configurationP, rndDir1, rndDir2, dirProbability); - dirProbability /= pSpecular; + //configurationP.energie /= (1-pSpecular) + //dirProbability /= (1 - pSpecular); + + //dirProbability = materialInfo->diffusion; + dirProbability = 1; + } else { nouvDirection = ReflectionLawsMLT::SolveSpecularPart(configurationP.direction, *materialInfo, faceNormal, configurationP, rndDir1, rndDir2, dirProbability); - dirProbability *= pSpecular; + //configurationP.energie /= pSpecular; + //dirProbability /= pSpecular; + + //dirProbability = (1-materialInfo->diffusion); + dirProbability = 1; + } //Calcul de la nouvelle direction de réflexion (en reprenant la célérité de propagation du son) configurationP.direction = nouvDirection * distanceToTravel; - configurationP.totalProbability *= dirProbability; + configurationP.reflectionOrder++; + //configurationP.energie *= dirProbability; + configurationP.probability.currentProbability *= dirProbability; return false; } @@ -301,8 +334,9 @@ void MLTCore::Follow(CONF_PARTICULE_MLT& propagationParticle) decimal distanceToTravel = propagationParticle.direction.length(); propagationParticle.ResetToInitialState(); - propagationParticle.totalProbability = 1; + propagationParticle.probability.currentProbability = 1; propagationParticle.totalDistance = 0; + propagationParticle.reflectionOrder = 0; for (int i = 0; i < propagationParticle.reflectionsNum; i++) { diff --git a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp index b61d82f6..94acebb9 100644 --- a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp +++ b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp @@ -307,7 +307,7 @@ void NextEventEstimationCore::FreeParticleTranslation(CONF_PARTICULE_AGH &config configurationP.position += translationVector; } -void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo,const vec3& faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace* faceInfo) +void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo,const vec3& faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace* faceInfo, std::vector* probability) { //Calculate and cast shadow rays for each (t_Recepteur_P* receiver in configurationTool->recepteur_p_List) @@ -352,8 +352,12 @@ void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t shadowRay.currentTetra = &sceneTetraMesh->tetraedres[receiver->indexTetra]; shadowRay.energie *= pow(densite_proba_absorption_atmospherique, timeStepNum); + if (probability != nullptr) { + //shadowRay.energie *= energy; + probability->push_back(energy); + } + shadowRays.push_back(shadowRay); - //if(probability != nullptr) *probability *= energy; } } } \ No newline at end of file diff --git a/src/spps-agh/calculation_cores/NextEventEstimationCore.h b/src/spps-agh/calculation_cores/NextEventEstimationCore.h index a81c0b56..576295e2 100644 --- a/src/spps-agh/calculation_cores/NextEventEstimationCore.h +++ b/src/spps-agh/calculation_cores/NextEventEstimationCore.h @@ -13,7 +13,7 @@ class NextEventEstimationCore : public CalculationCoreSPPS protected: virtual void Movement(CONF_PARTICULE_AGH &configurationP) override; virtual void FreeParticleTranslation(CONF_PARTICULE_AGH &configurationP, const vec3 &translationVector) override; - void GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq * materialInfo, const vec3 & faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace * faceInfo); + void GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq * materialInfo, const vec3 & faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace * faceInfo, std::vector* probability = nullptr); }; #endif \ No newline at end of file diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 9325c460..7c1658cf 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -63,12 +63,28 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH std::list shadowRays; double totalDistance = 0; - double totalProbability = 1; double weight = 0; int reflectionsNum = 0; int creationTime = 0; int largeStep = 0; + class + { + public: + double currentProbability = 1; + std::vector partialProbability; + + double totalProbability() + { + double totalProb = 0; + for(double val : partialProbability){ + totalProb += val; + } + return totalProb; + } + + } probability; + struct initialState { vec3 direction; diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index 19aeae80..21c82afe 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -86,7 +86,7 @@ class ReflectionLawsAGH : public ReflectionLaws * * Les vecteurs de réflexion retournés sont normalisés car une perte de précision se fait ressentir quand à la norme de ces vecteurs */ -class ReflectionLawsMLT : public ReflectionLaws +class ReflectionLawsMLT : public ReflectionLawsAGH { public: /** @@ -102,7 +102,7 @@ class ReflectionLawsMLT : public ReflectionLaws switch (materialInfo.reflectionLaw) { case REFLECTION_LAW_SPECULAR: - probability = (1 - materialInfo.diffusion); + probability = 1; return SpecularReflection(vectorDirection, faceNormal); case REFLECTION_LAW_LAMBERT: return BaseWnReflection(vectorDirection, faceNormal, 1, rnd1, rnd2, probability); @@ -128,7 +128,7 @@ class ReflectionLawsMLT : public ReflectionLaws case REFLECTION_LAW_PHONG: return PhongSpecularPart(vectorDirection, faceNormal, particuleInfo, materialInfo, rnd1, rnd2, probability); default: - probability = (1 - materialInfo.diffusion); + probability = 1; return SpecularReflection(vectorDirection, faceNormal); }; } @@ -138,7 +138,7 @@ class ReflectionLawsMLT : public ReflectionLaws static vec3 PhongSpecularPart(vec3 &vectorDirection, vec3 &faceNormal, CONF_PARTICULE_AGH& particuleInfo, t_Material_BFreq &material, double &rnd1, double &rnd2, double &probability) { - float n = powf(10, powf(-0.82662*material.diffusion, 3) + 1.5228); + float n = powf(10, -1.7234170470604733 * material.diffusion + 2.6245274102195886); vec3 specularDir = SpecularReflection(vectorDirection, faceNormal); Matrix3 rotationMatrix; @@ -146,12 +146,15 @@ class ReflectionLawsMLT : public ReflectionLaws vec3 target = BaseWnReflection(vectorDirection, faceNormal, n, rnd1, rnd2, probability); target = rotationMatrix*target; + target.normalize(); //test if reflected ray is pointing into ground if (target.dot(faceNormal) <= 0) { //if it is get new random direction - target = PhongSpecularPart(vectorDirection, faceNormal, particuleInfo, material, rnd1, rnd2, probability); + //if facing ground the ray is not possible to create, must fore ray reflection ???? + vec3 spec = SpecularReflection(target, faceNormal); + target = spec; } return target / target.length(); @@ -163,15 +166,10 @@ class ReflectionLawsMLT : public ReflectionLaws decimal theta = rnd1 * M_2PI; decimal phi = acos(pow((float)1 - rnd2, (float)(1. / (expo + 1.))));//pow((float)acos(1-GetRandValue()),(float)(1./(expo+1.))); - probability = ((expo + 2) / (2 * M_PI)) * powf(cos(phi),expo); + //probability = ((expo + 2) / (2 * M_PI)) * pow(cos(phi),expo); + //probability = (1/M_2PI) * pow(cos(phi), expo); + probability = pow(cos(phi), expo); return BaseUniformReflection(vecteurVitesse, faceNormal, theta, phi); } - - - static vec3 SpecularReflection(vec3 &vecteurVitesse, vec3 &faceNormal) - { - vec3 retVal = (vecteurVitesse - (faceNormal*(vecteurVitesse.dot(faceNormal)) * 2)); - return retVal / retVal.length(); - } }; \ No newline at end of file From 3956df484309dea57edf38665c9b7f224eb7fec4 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 25 Mar 2021 11:44:47 +0100 Subject: [PATCH 02/56] Initial source direction randomization [WIP] --- src/spps-agh/calculation_cores/MLTCore.cpp | 20 ++++++++++++++++++-- src/spps-agh/sppsNeeAGHTypes.h | 2 ++ src/spps/tools/dotdistribution.cpp | 15 ++++++++++++--- src/spps/tools/dotdistribution.h | 2 +- 4 files changed, 33 insertions(+), 6 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 382b662b..3de70b27 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -120,6 +120,14 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu bool MLTCore::RunInitialSeed(CONF_PARTICULE_MLT& inputParticle) { + // Recalculate source direction so it can be mutated + double rndDir1 = GetRandValue(); + double rndDir2 = GetRandValue(); + decimal distanceToTravel = inputParticle.direction.length(); + ParticleDistribution::GenSphereDistribution(inputParticle, distanceToTravel, &rndDir1, &rndDir2); + inputParticle.sourceDir1=rndDir1; + inputParticle.sourceDir2=rndDir2; + SetNextParticleCollision(inputParticle); //1er test de collision SetNextParticleCollisionWithObstructionElement(inputParticle); // Test de collision avec objet virtuel encombrant //Au premier pas de temps il faut enregistrer l'energie de la particule dans la maille courante @@ -179,8 +187,12 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML mutetedParticle.matAbsorbtions.clear(); mutetedParticle.largeStep = 1; mutetedParticle.weight = 0; + mutetedParticle.sourceDir1 = GetRandValue(); + mutetedParticle.sourceDir2 = GetRandValue(); }else{ + MutationGenerator(mutetedParticle.sourceDir1); + MutationGenerator(mutetedParticle.sourceDir2); for (int i = 0; i < input_particle.reflectionsNum; i++) { MutationGenerator(mutetedParticle.refDir1[i]); @@ -286,7 +298,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn //dirProbability /= (1 - pSpecular); //dirProbability = materialInfo->diffusion; - dirProbability = 1; + //dirProbability = 1; } else @@ -296,7 +308,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn //dirProbability /= pSpecular; //dirProbability = (1-materialInfo->diffusion); - dirProbability = 1; + //dirProbability = 1; } @@ -338,6 +350,10 @@ void MLTCore::Follow(CONF_PARTICULE_MLT& propagationParticle) propagationParticle.totalDistance = 0; propagationParticle.reflectionOrder = 0; + // set new source direction based on new random values + ParticleDistribution::GenSphereDistribution(propagationParticle, distanceToTravel, &propagationParticle.sourceDir1, &propagationParticle.sourceDir2); + + // follow reflections based on provided random values for (int i = 0; i < propagationParticle.reflectionsNum; i++) { double rndDir1 = propagationParticle.refDir1[i]; diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 7c1658cf..9d0cc315 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -60,6 +60,8 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH std::vector travelProbability; std::vector refDir1; std::vector refDir2; + double sourceDir1; + double sourceDir2; std::list shadowRays; double totalDistance = 0; diff --git a/src/spps/tools/dotdistribution.cpp b/src/spps/tools/dotdistribution.cpp index d2cc6c5a..cbe43328 100644 --- a/src/spps/tools/dotdistribution.cpp +++ b/src/spps/tools/dotdistribution.cpp @@ -3,15 +3,24 @@ -void ParticleDistribution::GenSphereDistribution(struct CONF_PARTICULE& particle,decimal vecNorm) +void ParticleDistribution::GenSphereDistribution(struct CONF_PARTICULE& particle,decimal vecNorm, double* rand1, double* rand2) { using namespace std; + if (rand1 == nullptr) { + double r1 = GetRandValue(); + rand1 = &r1; + } + if (rand2 == nullptr) { + double r2 = GetRandValue(); + rand2 = &r2; + } + // Choose a random longitude for the point - float phi = GetRandValue() * M_2PI; + float phi = *rand1 * M_2PI; // The z axis goes straight through the sphere and is chosen in random. - float z = GetRandValue() * 2.0f - 1.0f; + float z = *rand2 * 2.0f - 1.0f; // From that, we'll calculate the latitude this point will have on the // sphere's surface. diff --git a/src/spps/tools/dotdistribution.h b/src/spps/tools/dotdistribution.h index 5190786e..bc8a9333 100644 --- a/src/spps/tools/dotdistribution.h +++ b/src/spps/tools/dotdistribution.h @@ -18,7 +18,7 @@ class ParticleDistribution /** * Ajouter à la liste de particules */ - static void GenSphereDistribution(struct CONF_PARTICULE& particle,decimal vecNorm); + static void GenSphereDistribution(struct CONF_PARTICULE& particle,decimal vecNorm, double* rand1=nullptr, double* rand2= nullptr); static void GenXYDistribution(struct CONF_PARTICULE& particle,decimal vecNorm); static void GenXZDistribution(struct CONF_PARTICULE& particle,decimal vecNorm); static void GenYZDistribution(struct CONF_PARTICULE& particle,decimal vecNorm); From 21dceb2261ac22960cdb0f8987136448a1784c0c Mon Sep 17 00:00:00 2001 From: wbinek Date: Sun, 3 Oct 2021 14:24:01 +0200 Subject: [PATCH 03/56] Add warning about experimental code and its limitations --- src/isimpa/data_manager/projet.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/isimpa/data_manager/projet.cpp b/src/isimpa/data_manager/projet.cpp index 06cf1b96..88f0938e 100644 --- a/src/isimpa/data_manager/projet.cpp +++ b/src/isimpa/data_manager/projet.cpp @@ -694,6 +694,15 @@ void ProjectManager::RunCoreCalculation(Element* coreCalculation) return; } + if (coreCalculation->GetElementInfos().typeElement == Element::ELEMENT_TYPE_CORE_SPPSAGH) + { + E_Core_Core_Configuration* confCore = dynamic_cast(coreCalculation->GetElementByType(Element::ELEMENT_TYPE_CORE_CORE_CONFIG)); + int calcCode = confCore->GetListConfig("calculation_core"); + if (calcCode == 2) { + wxLogWarning(_("You are using very experimantal MLT calculation code. Many things are not working!\n - only omnidirectional source supported\n - only point receivers supported")); + } + } + if (coreCalculation->GetElementInfos().typeElement != Element::ELEMENT_TYPE_CORE_SPPSAGH) { std::vector lstMatRow, egroupeSurf; @@ -720,7 +729,6 @@ void ProjectManager::RunCoreCalculation(Element* coreCalculation) } } } - } //On active l'onglet des messages From 784d070a2885441830c6806542d0609c1d9464e8 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 14 Oct 2021 20:46:01 +0200 Subject: [PATCH 04/56] Add some comentes --- src/spps-agh/calculation_cores/MLTCore.cpp | 7 ++----- src/spps-agh/calculation_cores/MLTCore.h | 2 +- src/spps-agh/sppsNeeAGHTypes.h | 6 +++--- src/spps-agh/tools/dotreflectionAGH.h | 2 +- 4 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 3de70b27..bdddc908 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -65,7 +65,6 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) shadowRay = input_particle.shadowRays.front(); input_particle.shadowRays.pop_front(); //applicationTools.outputTool->NewParticule(confPart); - //shadowRay.energie *= input_particle.weight; shadowRay.energie *= input_particle.weight; Run(shadowRay); //applicationTools.outputTool->SaveParticule(); @@ -285,9 +284,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn //******************************* REFLECTION ***************************************// - // Choix de la méthode de reflexion en fonction de la valeur de diffusion - - double dirProbability; + double dirProbability; // probability of sending ray in given direction (what about scattering coefficient)?? //Get direction for diffuse or specular part based on material info pSpecular = 1; @@ -316,7 +313,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn configurationP.direction = nouvDirection * distanceToTravel; configurationP.reflectionOrder++; //configurationP.energie *= dirProbability; - configurationP.probability.currentProbability *= dirProbability; + configurationP.probability.currentProbability *= dirProbability; //optionally multiply by scattering coefficient (althouth current implementation forces more scattered rays?)??? return false; } diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 836b32f6..33dab48c 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -19,7 +19,7 @@ class MLTCore : public NextEventEstimationCore int mutation_number = 25; protected: - std::list Seeds; + std::list Seeds; void Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_MLT& mutetedParticle) const; bool MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rndDir1, double rndDir2, double rndTravelProbability, double rndMatAbsorbtion, float deltaT, float distanceToTravel); diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 9d0cc315..d97460d0 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -73,10 +73,10 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH class { public: - double currentProbability = 1; - std::vector partialProbability; + double currentProbability = 1; //probability of generating ray path up to some point + std::vector partialProbability; //probability of generating shadow ray - double totalProbability() + double totalProbability() //sum of shadow rays probability - total ray contribution { double totalProb = 0; for(double val : partialProbability){ diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index 21c82afe..51f1e6cb 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -168,7 +168,7 @@ class ReflectionLawsMLT : public ReflectionLawsAGH //probability = ((expo + 2) / (2 * M_PI)) * pow(cos(phi),expo); //probability = (1/M_2PI) * pow(cos(phi), expo); - probability = pow(cos(phi), expo); + probability = pow(cos(phi), expo); //normalization factor excluded, to aviod large numbers return BaseUniformReflection(vecteurVitesse, faceNormal, theta, phi); } From 23386542661fb2b7157b5af23171f274bee85ae1 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 14 Oct 2021 20:46:25 +0200 Subject: [PATCH 05/56] Add new MLT parameters --- .../tree_core/e_core_spps_aghcore_advanced.h | 15 ++++++++ .../data_manager/core_configurationAGH.cpp | 36 +++++++++++++++++++ .../data_manager/core_configurationAGH.h | 11 ++++++ 3 files changed, 62 insertions(+) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index dffb6f41..80bd6867 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -223,6 +223,9 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element if (!this->IsPropertyExist("mutationNumber")) { InitBaseMLTProperties(this); } + if (!this->IsPropertyExist("useRejectionSampling")) { + InitAdditionalMLTProperties(this); + } } E_Core_SppsNee_AGH_advanced_MLT(Element* parent) @@ -258,6 +261,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element { InitRandomSeed(this); InitBaseMLTProperties(this); + InitAdditionalMLTProperties(this); } void InitRandomSeed(Element* confCore) { @@ -269,6 +273,17 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Mutation number"), 25, true, false, true, 0, 0); } + void InitAdditionalMLTProperties(Element* confCore) { + confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use rejection sampling"), true, true); + confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Seed to sample ratio"), 0.5, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use feature scaling"), true, true); + confCore->AppendPropertyBool("mutateDirectSound", wxTRANSLATE("Mutate direct sound"), true, true); + confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); + confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), true, true); + confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); + confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), true, true); + confCore->AppendPropertyBool("costAirAbs", wxTRANSLATE("Cost - use air absorption"), true, true); + } }; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 9090959e..b17f9db6 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -198,6 +198,42 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, advancedNode->GetProperty("specularReflProbabilityMulti").ToFloat()); SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, advancedNode->GetProperty("mutationNumber").ToInt()); SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, advancedNode->GetProperty("largeStepProb").ToFloat()); + + //Load new MLT parameters, do check on them, write to console if not implemented + if (advancedNode->GetProperty("useRejectionSampling").ToInt()) { + std::cout << "Rejection sampling not implemented yet!!" << std::endl; + } + if (advancedNode->GetProperty("useFeatureScaling").ToInt()) { + std::cout << "Feature scaling not implemented yet!!" << std::endl; + } + if (advancedNode->GetProperty("mutateDirectSound").ToInt()) { + std::cout << "Direct sound mutation - currently enabled only, limited functionality, omni source only!!" << std::endl; + } + if (advancedNode->GetProperty("costMatAbsorption").ToInt()) { + std::cout << "Cost - absorption - currently permanently disabled" << std::endl; + } + if (advancedNode->GetProperty("costMatDiffusion").ToInt()) { + std::cout << "Cost - diffusion - currently permanently disabled" << std::endl; + } + if (advancedNode->GetProperty("costShadowRays").ToInt()) { + std::cout << "Cost - SR - currently permanently enabled" << std::endl; + } + if (advancedNode->GetProperty("costReflDir").ToInt()) { + std::cout << "Cost - reflection dir. - currently permanently enabled" << std::endl; + } + if (advancedNode->GetProperty("costAirAbs").ToInt()) { + std::cout << "Cost - air absoprtion - currently permanently disabled" << std::endl; + } + + SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, advancedNode->GetProperty("useRejectionSampling").ToInt()); + SetConfigInformation(IPROP_MLT_DO_FEATURE_SCALING, advancedNode->GetProperty("useFeatureScaling").ToInt()); + SetConfigInformation(IPROP_MLT_DO_MUTATE_DIRECT_SOUND, advancedNode->GetProperty("mutateDirectSound").ToInt()); + SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, advancedNode->GetProperty("costMatAbsorption").ToInt()); + SetConfigInformation(IPROP_MLT_COST_MAT_DIFFUSION, advancedNode->GetProperty("costMatDiffusion").ToInt()); + SetConfigInformation(IPROP_MLT_COST_SHADOW_RAYS, advancedNode->GetProperty("costShadowRays").ToInt()); + SetConfigInformation(IPROP_MLT_COST_REFL_DIR, advancedNode->GetProperty("costReflDir").ToInt()); + SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, advancedNode->GetProperty("costAirAbs").ToInt()); + SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, advancedNode->GetProperty("seedSampleRatio").ToFloat()); } } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index 288e3168..d5837412 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -30,6 +30,7 @@ class Core_ConfigurationAGH : public Core_Configuration FPROP_MLT_SPECULAR_REFL_PROB, FPROP_MLT_MUTATION_NUMBER, FPROP_NEE_SHADOWRAY_PROB, + FPROP_MLT_SEED_SAMPLE_RATIO, }; enum NEW_IPROP @@ -40,6 +41,16 @@ class Core_ConfigurationAGH : public Core_Configuration IPROP_MAP_MIN_REFL, /*!< Min reflection order for map calculation*/ IPROP_SKIP_DIRECT_SOUND_CALC, IPROP_MAP_MAX_REFL, /*!< Max reflection order for map calculation*/ + + //MLT related boolean properties + IPROP_MLT_DO_REJECTION_SAMPLING, + IPROP_MLT_DO_FEATURE_SCALING, + IPROP_MLT_DO_MUTATE_DIRECT_SOUND, + IPROP_MLT_COST_MAT_ABSORPTION, + IPROP_MLT_COST_MAT_DIFFUSION, + IPROP_MLT_COST_SHADOW_RAYS, + IPROP_MLT_COST_REFL_DIR, + IPROP_MLT_COST_AIR_ABSORPTION }; /** * Initialisation des paramètres du coeur de calcul à partir d'un fichier XML From 811e2ad69785b53ce0ec576c744af4a6b7637748 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 15 Oct 2021 12:44:58 +0200 Subject: [PATCH 06/56] Implement most of MLT parameters as variables --- .../tree_core/e_core_spps_aghcore_advanced.h | 5 +- src/spps-agh/calculation_cores/MLTCore.cpp | 44 +++++++----- src/spps-agh/calculation_cores/MLTCore.h | 8 +++ .../data_manager/core_configurationAGH.cpp | 68 +++++++++++++------ .../data_manager/core_configurationAGH.h | 3 +- 5 files changed, 87 insertions(+), 41 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 80bd6867..f732e477 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -276,13 +276,14 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element void InitAdditionalMLTProperties(Element* confCore) { confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use rejection sampling"), true, true); confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Seed to sample ratio"), 0.5, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use feature scaling"), true, true); confCore->AppendPropertyBool("mutateDirectSound", wxTRANSLATE("Mutate direct sound"), true, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); - confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), true, true); + confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), true, true); - confCore->AppendPropertyBool("costAirAbs", wxTRANSLATE("Cost - use air absorption"), true, true); + confCore->AppendPropertyBool("costAirAbs", wxTRANSLATE("Cost - use air absorption"), false, true); } }; diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index bdddc908..a10688fe 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -14,6 +14,14 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT pLarge = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_LARGE_STEP_PROB); pSpecular = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SPECULAR_REFL_PROB); mutation_number = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); + reflectionPenaltyMultip = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_REFL_PENALTY_MULTIP); + + costMatAbsorption = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_ABSORPTION); + costMatDiffusion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_DIFFUSION); + costReflDir = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_REFL_DIR); + costShadowRays = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_SHADOW_RAYS); + costAirAbsoprtion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_AIR_ABSORPTION); + mutateSource = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_MUTATE_SOURCE); }; bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) @@ -224,7 +232,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn // return true; //} configurationP.energie *= travelProbability; - //configurationP.probability.currentProbability *= travelProbability; + if(costAirAbsoprtion) configurationP.probability.currentProbability *= travelProbability; configurationP.totalDistance += distance; @@ -249,7 +257,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn //} configurationP.energie *= (1 - materialInfo->absorption); - //configurationP.probability.currentProbability *= (1 - materialInfo->absorption); + if (costMatAbsorption) configurationP.probability.currentProbability *= (1 - materialInfo->absorption); //******************************* UPDATE TIME STEP NUMBER ***********************************// double timestepNum = configurationP.totalDistance / configurationP.direction.length(); @@ -277,8 +285,14 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn if (configurationP.reflectionOrder > 0) { - for (double srProb : srProbabilities) { - configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability * srProb); + if(costShadowRays) // If shadow ray cost is enabled adds the currentProb * srProb for each Shadow Ray + { + for (double srProb : srProbabilities) { + configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability * srProb); + } + } + else { // If shadow ray cost is disablerd adds the currentProb once (this option is rather useless) + configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability); } } @@ -291,29 +305,25 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn if (materialInfo->diffusion == 1 || pSpecular*((1 - rndMatAbsorbtion)) < materialInfo->diffusion) { nouvDirection = ReflectionLawsMLT::SolveDiffusePart(configurationP.direction, *materialInfo, faceNormal, configurationP, rndDir1, rndDir2, dirProbability); + if (costMatDiffusion) configurationP.probability.currentProbability *= materialInfo->diffusion; //configurationP.energie /= (1-pSpecular) - //dirProbability /= (1 - pSpecular); - - //dirProbability = materialInfo->diffusion; - //dirProbability = 1; - + //dirProbability /= (1 - pSpecular); } else { nouvDirection = ReflectionLawsMLT::SolveSpecularPart(configurationP.direction, *materialInfo, faceNormal, configurationP, rndDir1, rndDir2, dirProbability); + if (costMatDiffusion) configurationP.probability.currentProbability *= (1 - materialInfo->diffusion); //configurationP.energie /= pSpecular; - //dirProbability /= pSpecular; - - //dirProbability = (1-materialInfo->diffusion); - //dirProbability = 1; - + //dirProbability /= pSpecular; } //Calcul de la nouvelle direction de réflexion (en reprenant la célérité de propagation du son) configurationP.direction = nouvDirection * distanceToTravel; configurationP.reflectionOrder++; //configurationP.energie *= dirProbability; - configurationP.probability.currentProbability *= dirProbability; //optionally multiply by scattering coefficient (althouth current implementation forces more scattered rays?)??? + if(costReflDir) configurationP.probability.currentProbability *= dirProbability; + + configurationP.probability.currentProbability *= reflectionPenaltyMultip; return false; } @@ -348,7 +358,9 @@ void MLTCore::Follow(CONF_PARTICULE_MLT& propagationParticle) propagationParticle.reflectionOrder = 0; // set new source direction based on new random values - ParticleDistribution::GenSphereDistribution(propagationParticle, distanceToTravel, &propagationParticle.sourceDir1, &propagationParticle.sourceDir2); + if (mutateSource) { + ParticleDistribution::GenSphereDistribution(propagationParticle, distanceToTravel, &propagationParticle.sourceDir1, &propagationParticle.sourceDir2); + } // follow reflections based on provided random values for (int i = 0; i < propagationParticle.reflectionsNum; i++) diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 33dab48c..8d46f416 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -17,6 +17,14 @@ class MLTCore : public NextEventEstimationCore float pLarge = 0.5; float pSpecular = 0.5; int mutation_number = 25; + float reflectionPenaltyMultip = 0.95; + + bool costMatAbsorption = true; + bool costMatDiffusion = false; + bool costReflDir = true; + bool costShadowRays = true; + bool costAirAbsoprtion = false; + bool mutateSource = true; protected: std::list Seeds; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index b17f9db6..1524becb 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -195,44 +195,68 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Random seed has been set; then multi-thread has been desactivated." << std::endl; } SetConfigInformation(I_PROP_RANDOM_SEED, seed); - SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, advancedNode->GetProperty("specularReflProbabilityMulti").ToFloat()); - SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, advancedNode->GetProperty("mutationNumber").ToInt()); - SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, advancedNode->GetProperty("largeStepProb").ToFloat()); - //Load new MLT parameters, do check on them, write to console if not implemented + + //Load MLT parameters, do check on them, write to console if not implemented + int mutationNumber = advancedNode->GetProperty("mutationNumber").ToInt(); + float largeStepProb = advancedNode->GetProperty("largeStepProb").ToFloat(); + float specularReflProbabilityMulti = advancedNode->GetProperty("specularReflProbabilityMulti").ToFloat(); + + float reflPenaltyMultip = advancedNode->GetProperty("reflectionPenalty").ToFloat(); + int costMatAbsorption = advancedNode->GetProperty("costMatAbsorption").ToInt(); + int costMatDiffusion = advancedNode->GetProperty("costMatDiffusion").ToInt(); + int costReflDir = advancedNode->GetProperty("costReflDir").ToInt(); + int costShadowRays = advancedNode->GetProperty("costShadowRays").ToInt(); + int costAirAbsoprtion = advancedNode->GetProperty("costAirAbs").ToInt(); + int mutateSource = advancedNode->GetProperty("mutateSource").ToInt(); + + + std::cout << "Mutation number: " << mutationNumber << std::endl; + std::cout << "Large step probability: " << largeStepProb << std::endl << + "Currently large steps are disabled/not implemented" << std::endl; + std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << + "Currently spec. ref. multiplier is disabled" << std::endl; + std::cout << "Refection penalty multiplier " << reflPenaltyMultip << std::endl; + + if (advancedNode->GetProperty("useRejectionSampling").ToInt()) { std::cout << "Rejection sampling not implemented yet!!" << std::endl; } if (advancedNode->GetProperty("useFeatureScaling").ToInt()) { std::cout << "Feature scaling not implemented yet!!" << std::endl; } - if (advancedNode->GetProperty("mutateDirectSound").ToInt()) { - std::cout << "Direct sound mutation - currently enabled only, limited functionality, omni source only!!" << std::endl; + if (mutateSource) { + std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; } - if (advancedNode->GetProperty("costMatAbsorption").ToInt()) { - std::cout << "Cost - absorption - currently permanently disabled" << std::endl; + if (costMatAbsorption) { + std::cout << "Cost - absorption: "<< costMatAbsorption << std::endl; } - if (advancedNode->GetProperty("costMatDiffusion").ToInt()) { - std::cout << "Cost - diffusion - currently permanently disabled" << std::endl; + if (costMatDiffusion) { + std::cout << "Cost - diffusion: " << costMatDiffusion << std::endl; } - if (advancedNode->GetProperty("costShadowRays").ToInt()) { - std::cout << "Cost - SR - currently permanently enabled" << std::endl; + if (costShadowRays) { + std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; } - if (advancedNode->GetProperty("costReflDir").ToInt()) { - std::cout << "Cost - reflection dir. - currently permanently enabled" << std::endl; + if (costReflDir) { + std::cout << "Cost - reflection dir.: "<< costReflDir << std::endl; } - if (advancedNode->GetProperty("costAirAbs").ToInt()) { - std::cout << "Cost - air absoprtion - currently permanently disabled" << std::endl; + if (costAirAbsoprtion) { + std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; } + + SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); + SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, largeStepProb); + SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, specularReflProbabilityMulti); + SetConfigInformation(FPROP_MLT_REFL_PENALTY_MULTIP, reflPenaltyMultip); SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, advancedNode->GetProperty("useRejectionSampling").ToInt()); SetConfigInformation(IPROP_MLT_DO_FEATURE_SCALING, advancedNode->GetProperty("useFeatureScaling").ToInt()); - SetConfigInformation(IPROP_MLT_DO_MUTATE_DIRECT_SOUND, advancedNode->GetProperty("mutateDirectSound").ToInt()); - SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, advancedNode->GetProperty("costMatAbsorption").ToInt()); - SetConfigInformation(IPROP_MLT_COST_MAT_DIFFUSION, advancedNode->GetProperty("costMatDiffusion").ToInt()); - SetConfigInformation(IPROP_MLT_COST_SHADOW_RAYS, advancedNode->GetProperty("costShadowRays").ToInt()); - SetConfigInformation(IPROP_MLT_COST_REFL_DIR, advancedNode->GetProperty("costReflDir").ToInt()); - SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, advancedNode->GetProperty("costAirAbs").ToInt()); + SetConfigInformation(IPROP_MLT_DO_MUTATE_SOURCE, mutateSource); + SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, costMatAbsorption); + SetConfigInformation(IPROP_MLT_COST_MAT_DIFFUSION, costMatDiffusion); + SetConfigInformation(IPROP_MLT_COST_SHADOW_RAYS, costShadowRays); + SetConfigInformation(IPROP_MLT_COST_REFL_DIR, costReflDir); + SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, costAirAbsoprtion); SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, advancedNode->GetProperty("seedSampleRatio").ToFloat()); } } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index d5837412..b6509dd4 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -31,6 +31,7 @@ class Core_ConfigurationAGH : public Core_Configuration FPROP_MLT_MUTATION_NUMBER, FPROP_NEE_SHADOWRAY_PROB, FPROP_MLT_SEED_SAMPLE_RATIO, + FPROP_MLT_REFL_PENALTY_MULTIP, }; enum NEW_IPROP @@ -45,7 +46,7 @@ class Core_ConfigurationAGH : public Core_Configuration //MLT related boolean properties IPROP_MLT_DO_REJECTION_SAMPLING, IPROP_MLT_DO_FEATURE_SCALING, - IPROP_MLT_DO_MUTATE_DIRECT_SOUND, + IPROP_MLT_DO_MUTATE_SOURCE, IPROP_MLT_COST_MAT_ABSORPTION, IPROP_MLT_COST_MAT_DIFFUSION, IPROP_MLT_COST_SHADOW_RAYS, From 2f17eae3fcea0114830f7072db8195eab995f39e Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 15 Oct 2021 20:03:08 +0200 Subject: [PATCH 07/56] Rename MLT parameter --- .../data_manager/tree_core/e_core_spps_aghcore_advanced.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index f732e477..cfaedf2c 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -278,7 +278,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Seed to sample ratio"), 0.5, false, 2, true, true, 1, 0, true); confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use feature scaling"), true, true); - confCore->AppendPropertyBool("mutateDirectSound", wxTRANSLATE("Mutate direct sound"), true, true); + confCore->AppendPropertyBool("mutateSource", wxTRANSLATE("Mutate source"), true, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); From ab647193b494b150ad4f1b1955e10ae082cc43d9 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 15 Oct 2021 20:04:56 +0200 Subject: [PATCH 08/56] Add rejection sampling and seed choice, fix progress for MLT --- .../input_output/progressionInfo.h | 6 ++ src/spps-agh/calculation_cores/MLTCore.cpp | 82 +++++++++++++++---- src/spps-agh/calculation_cores/MLTCore.h | 7 ++ .../data_manager/core_configurationAGH.cpp | 34 +++----- src/spps-agh/sppsNeeAGH.cpp | 10 ++- 5 files changed, 101 insertions(+), 38 deletions(-) diff --git a/src/lib_interface/input_output/progressionInfo.h b/src/lib_interface/input_output/progressionInfo.h index 4129a2bf..9cce3f05 100644 --- a/src/lib_interface/input_output/progressionInfo.h +++ b/src/lib_interface/input_output/progressionInfo.h @@ -73,6 +73,12 @@ class progressOperation { PushProgression(1.f); } + + void Push(float increment) + { + PushProgression(increment); + } + protected: //friend class progressOperation; diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index a10688fe..69c93b5e 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -15,7 +15,9 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT pSpecular = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SPECULAR_REFL_PROB); mutation_number = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); reflectionPenaltyMultip = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_REFL_PENALTY_MULTIP); - + seedSampleRatio = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); + + rejectionSampling = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_REJECTION_SAMPLING); costMatAbsorption = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_ABSORPTION); costMatDiffusion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_DIFFUSION); costReflDir = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_REFL_DIR); @@ -33,16 +35,15 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) { if (!configurationP.isShadowRay) { - CONF_PARTICULE_MLT inputParticle(configurationP); inputParticle.SetInitialState(inputParticle); inputParticle.weight = 0; - int totalParticleNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); + int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); RunInitialSeed(inputParticle); - Seeds.push_back(inputParticle); - b += inputParticle.probability.totalProbability() / totalParticleNumber; + InitialSeeds.push_back(inputParticle); + b += inputParticle.probability.totalProbability() / totalSeedNumber; } else { @@ -52,6 +53,52 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) return true; } +void MLTCore::ChooseSeeds() { + InitialSeedsSize = InitialSeeds.size(); + + // If seed chosing disabled + if (seedSampleRatio == 0) { + Seeds = InitialSeeds; + std::cout << "Seed choosing disabled, copied all initial seeds" << std::endl; + return; + } + + // Seed choosing algorithm + int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); + float N_target = seedSampleRatio * totalSeedNumber; + + int potentialRepeatedRays = 0; + int totalNumberOfRepeats = 0; + + for (CONF_PARTICULE_MLT seed : InitialSeeds) { + float currentScore = seed.probability.totalProbability() / b * (N_target/totalSeedNumber); + if (currentScore > 1) { + potentialRepeatedRays++; + } + + while (currentScore > 1) { + totalNumberOfRepeats++; + Seeds.push_back(seed); + currentScore -= 1; + } + + float rndVal = GetRandValue(); + if (currentScore > rndVal) { + Seeds.push_back(seed); + } + } + SeedsSize = Seeds.size(); + + std::cout << "Final seed count: " << SeedsSize << std::endl; + std::cout << "Repeated rays: " << potentialRepeatedRays << std::endl; //both values are slightly wrong + std::cout << "Total repeats: " << totalNumberOfRepeats << std::endl; //both values are slightly wrong + if ((float)totalNumberOfRepeats / SeedsSize > 0.5) { + std::cout << "WARNING Repeated rays are " << ((float)totalNumberOfRepeats / SeedsSize) << + " of the total sample count!!" << std::endl << + "Consider lower sample to seed ratio" << std::endl; + } +} + bool MLTCore::SeedIsEmpty() { return Seeds.empty(); @@ -74,6 +121,7 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) input_particle.shadowRays.pop_front(); //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; + if (!seedSampleRatio == 0) shadowRay.energie *= InitialSeedsSize / SeedsSize; // Compensate energy for seed selection Run(shadowRay); //applicationTools.outputTool->SaveParticule(); } @@ -97,28 +145,32 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); + + if (rejectionSampling) { + // Rejection sampling + inputParticle.weight += b / mutation_number / inputParticle.probability.totalProbability() * (1 - ratio); + mutatedParticle.weight += b / mutation_number / mutatedParticle.probability.totalProbability() * ratio; + } + else { + // No rejection sampling + inputParticle.weight = b / inputParticle.probability.totalProbability(); + mutatedParticle.weight = b / mutatedParticle.probability.totalProbability(); + } - // Rejection sampling - //inputParticle.weight += b / mutation_number / inputParticle.probability.totalProbability() * (1- ratio); - //mutatedParticle.weight += b / mutation_number / mutatedParticle.probability.totalProbability() * ratio; - - // No rejection sampling - inputParticle.weight = b / inputParticle.probability.totalProbability(); - mutatedParticle.weight = b / mutatedParticle.probability.totalProbability(); if(mutatedParticle.largeStep==1) { - //CastShadowRays(inputParticle); + if (rejectionSampling) CastShadowRays(inputParticle); inputParticle = mutatedParticle; inputParticle.largeStep = 0; } else if(ratio > GetRandValue()) { - //CastShadowRays(inputParticle); + if (rejectionSampling) CastShadowRays(inputParticle); inputParticle = mutatedParticle; }else { - //CastShadowRays(mutatedParticle); + if (rejectionSampling) CastShadowRays(mutatedParticle); } mutation_counter++; } diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 8d46f416..80350e19 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -11,6 +11,7 @@ class MLTCore : public NextEventEstimationCore bool Run(CONF_PARTICULE_AGH configurationP) override; bool SeedIsEmpty(); void RunMutation(); + void ChooseSeeds(); float b = 0; //bias value - auto calculated - start with 0!; @@ -18,7 +19,9 @@ class MLTCore : public NextEventEstimationCore float pSpecular = 0.5; int mutation_number = 25; float reflectionPenaltyMultip = 0.95; + float seedSampleRatio = 0; + bool rejectionSampling = true; bool costMatAbsorption = true; bool costMatDiffusion = false; bool costReflDir = true; @@ -27,8 +30,12 @@ class MLTCore : public NextEventEstimationCore bool mutateSource = true; protected: + std::list InitialSeeds; std::list Seeds; + int InitialSeedsSize; + int SeedsSize; + void Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_MLT& mutetedParticle) const; bool MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rndDir1, double rndDir2, double rndTravelProbability, double rndMatAbsorbtion, float deltaT, float distanceToTravel); void CastShadowRays(CONF_PARTICULE_MLT& input_particle); diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 1524becb..a1552e9b 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -209,14 +209,24 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int costShadowRays = advancedNode->GetProperty("costShadowRays").ToInt(); int costAirAbsoprtion = advancedNode->GetProperty("costAirAbs").ToInt(); int mutateSource = advancedNode->GetProperty("mutateSource").ToInt(); + float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); - + largeStepProb = 0; // to be deleted after proper large step implementation, currently fixed to avoid errors std::cout << "Mutation number: " << mutationNumber << std::endl; std::cout << "Large step probability: " << largeStepProb << std::endl << "Currently large steps are disabled/not implemented" << std::endl; std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << "Currently spec. ref. multiplier is disabled" << std::endl; - std::cout << "Refection penalty multiplier " << reflPenaltyMultip << std::endl; + std::cout << "Refection penalty multiplier :" << reflPenaltyMultip << std::endl; + std::cout << "Seed to sample rato :" << seedSampleRatio << std::endl; + if(seedSampleRatio==0) std::cout << "Seed selection disabled!" << std::endl; + + std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; + std::cout << "Cost - absorption: " << costMatAbsorption << std::endl; + std::cout << "Cost - diffusion: " << costMatDiffusion << std::endl; + std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; + std::cout << "Cost - reflection dir.: " << costReflDir << std::endl; + std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; if (advancedNode->GetProperty("useRejectionSampling").ToInt()) { @@ -225,24 +235,6 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { if (advancedNode->GetProperty("useFeatureScaling").ToInt()) { std::cout << "Feature scaling not implemented yet!!" << std::endl; } - if (mutateSource) { - std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; - } - if (costMatAbsorption) { - std::cout << "Cost - absorption: "<< costMatAbsorption << std::endl; - } - if (costMatDiffusion) { - std::cout << "Cost - diffusion: " << costMatDiffusion << std::endl; - } - if (costShadowRays) { - std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; - } - if (costReflDir) { - std::cout << "Cost - reflection dir.: "<< costReflDir << std::endl; - } - if (costAirAbsoprtion) { - std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; - } SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, largeStepProb); @@ -257,7 +249,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_COST_SHADOW_RAYS, costShadowRays); SetConfigInformation(IPROP_MLT_COST_REFL_DIR, costReflDir); SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, costAirAbsoprtion); - SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, advancedNode->GetProperty("seedSampleRatio").ToFloat()); + SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, seedSampleRatio); } } diff --git a/src/spps-agh/sppsNeeAGH.cpp b/src/spps-agh/sppsNeeAGH.cpp index 616c558d..ab2d12fc 100644 --- a/src/spps-agh/sppsNeeAGH.cpp +++ b/src/spps-agh/sppsNeeAGH.cpp @@ -89,7 +89,11 @@ void runSourceCalculation( progressOperation* parentOperation, t_ToolBox& applic } int progresSteps = quandparticules; - if (*applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::I_PROP_CALCULATION_CORE_SELLECTION) == MLT) progresSteps *= 2; + if (*applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::I_PROP_CALCULATION_CORE_SELLECTION) == MLT) { + int mutationNumber = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); + float seedSampleRatio = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); + progresSteps += seedSampleRatio*quandparticules* mutationNumber; + } progressOperation thisSrcOperation(parentOperation, progresSteps); @@ -171,6 +175,8 @@ void runSourceCalculation( progressOperation* parentOperation, t_ToolBox& applic if(*applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::I_PROP_CALCULATION_CORE_SELLECTION)==MLT) { + int mutationNumber = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); + static_cast(applicationTools.calculationTool)->ChooseSeeds(); // choose seeds from the initial seed pool while (!static_cast(applicationTools.calculationTool)->SeedIsEmpty()) { static_cast(applicationTools.calculationTool)->RunMutation(); @@ -179,7 +185,7 @@ void runSourceCalculation( progressOperation* parentOperation, t_ToolBox& applic boost::mutex::scoped_lock lock(amutex); #endif applicationTools.mainProgressionOutput->OutputCurrentProgression(); - thisSrcOperation.Next(); + thisSrcOperation.Push(mutationNumber); } } } From 1d1ddbfd5c89baad78c16a2d23e0e7d8394a8a0e Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 12:42:43 +0200 Subject: [PATCH 09/56] Fix very rare nan in BaseWnReflection probability calculation --- src/spps-agh/tools/dotreflectionAGH.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index 51f1e6cb..db818579 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -1,7 +1,6 @@ //#include "coreTypes.h" #include "tools/dotreflection.h" - /** * @brief Cette classe regroupe les méthodes de réfléxion de particules sur les surfaces * @@ -169,6 +168,9 @@ class ReflectionLawsMLT : public ReflectionLawsAGH //probability = ((expo + 2) / (2 * M_PI)) * pow(cos(phi),expo); //probability = (1/M_2PI) * pow(cos(phi), expo); probability = pow(cos(phi), expo); //normalization factor excluded, to aviod large numbers + if (isnan(probability)) [[unlikely]] { + probability = 0; + } return BaseUniformReflection(vecteurVitesse, faceNormal, theta, phi); } From ff94114ecc723c98a2bc37a33c9b77f2d6611942 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 13:41:17 +0200 Subject: [PATCH 10/56] Fix rejection sampling parameter loading --- src/spps-agh/data_manager/core_configurationAGH.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index a1552e9b..038bcdab 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -210,6 +210,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int costAirAbsoprtion = advancedNode->GetProperty("costAirAbs").ToInt(); int mutateSource = advancedNode->GetProperty("mutateSource").ToInt(); float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); + int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); largeStepProb = 0; // to be deleted after proper large step implementation, currently fixed to avoid errors std::cout << "Mutation number: " << mutationNumber << std::endl; @@ -222,16 +223,13 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { if(seedSampleRatio==0) std::cout << "Seed selection disabled!" << std::endl; std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; + std::cout << "Rejecton sampling: " << rejectionSampling << std::endl; std::cout << "Cost - absorption: " << costMatAbsorption << std::endl; std::cout << "Cost - diffusion: " << costMatDiffusion << std::endl; std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; std::cout << "Cost - reflection dir.: " << costReflDir << std::endl; std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; - - if (advancedNode->GetProperty("useRejectionSampling").ToInt()) { - std::cout << "Rejection sampling not implemented yet!!" << std::endl; - } if (advancedNode->GetProperty("useFeatureScaling").ToInt()) { std::cout << "Feature scaling not implemented yet!!" << std::endl; } @@ -241,7 +239,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, specularReflProbabilityMulti); SetConfigInformation(FPROP_MLT_REFL_PENALTY_MULTIP, reflPenaltyMultip); - SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, advancedNode->GetProperty("useRejectionSampling").ToInt()); + SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, rejectionSampling); SetConfigInformation(IPROP_MLT_DO_FEATURE_SCALING, advancedNode->GetProperty("useFeatureScaling").ToInt()); SetConfigInformation(IPROP_MLT_DO_MUTATE_SOURCE, mutateSource); SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, costMatAbsorption); From 7cc79890007647f8bc1c5d8e4fe54066d47eb625 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 16:59:34 +0200 Subject: [PATCH 11/56] Rename some MLT parameters --- .../tree_core/e_core_spps_aghcore_advanced.h | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index cfaedf2c..c050494b 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -269,16 +269,18 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element } void InitBaseMLTProperties(Element* confCore) { - confCore->AppendPropertyDecimal("largeStepProb", wxTRANSLATE("Large step probability"), 0.5, false, 2, true, true, 1, 0, true); - confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); - confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Mutation number"), 25, true, false, true, 0, 0); + confCore->AppendPropertyDecimal("largeStepProb", wxTRANSLATE("Tune - Large step probability"), 0.5, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Tune - Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); + confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Sampling - Mutation number"), 25, true, false, true, 0, 0); } void InitAdditionalMLTProperties(Element* confCore) { - confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use rejection sampling"), true, true); - confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Seed to sample ratio"), 0.5, false, 2, true, true, 1, 0, true); - confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); - confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use feature scaling"), true, true); - confCore->AppendPropertyBool("mutateSource", wxTRANSLATE("Mutate source"), true, true); + confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use - Rejection sampling"), true, true); + confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use - Feature scaling"), true, true); + confCore->AppendPropertyBool("mutateSource", wxTRANSLATE("Use - Mutate source"), true, true); + + confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); + confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); From 4eece4daae7a996010c71d0875c6ac3ce4bc85d3 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 16:59:58 +0200 Subject: [PATCH 12/56] [WIP] Large step sampling --- src/spps-agh/calculation_cores/MLTCore.cpp | 44 ++++++++++++------- .../data_manager/core_configurationAGH.cpp | 1 - 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 69c93b5e..4bd725f7 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -67,17 +67,18 @@ void MLTCore::ChooseSeeds() { int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); float N_target = seedSampleRatio * totalSeedNumber; - int potentialRepeatedRays = 0; - int totalNumberOfRepeats = 0; + int individualRays = 0; + bool accounted = false; for (CONF_PARTICULE_MLT seed : InitialSeeds) { + accounted = false; float currentScore = seed.probability.totalProbability() / b * (N_target/totalSeedNumber); if (currentScore > 1) { - potentialRepeatedRays++; + individualRays++; + accounted = true; } while (currentScore > 1) { - totalNumberOfRepeats++; Seeds.push_back(seed); currentScore -= 1; } @@ -85,16 +86,17 @@ void MLTCore::ChooseSeeds() { float rndVal = GetRandValue(); if (currentScore > rndVal) { Seeds.push_back(seed); + if (!accounted) individualRays++; } } SeedsSize = Seeds.size(); + float avgRayRepeat = (float)SeedsSize / individualRays; std::cout << "Final seed count: " << SeedsSize << std::endl; - std::cout << "Repeated rays: " << potentialRepeatedRays << std::endl; //both values are slightly wrong - std::cout << "Total repeats: " << totalNumberOfRepeats << std::endl; //both values are slightly wrong - if ((float)totalNumberOfRepeats / SeedsSize > 0.5) { - std::cout << "WARNING Repeated rays are " << ((float)totalNumberOfRepeats / SeedsSize) << - " of the total sample count!!" << std::endl << + std::cout << "Individual rays: " << individualRays << std::endl; //both values are slightly wrong + std::cout << "Average ray repeat: " << avgRayRepeat << std::endl; //both values are slightly wrong + if (avgRayRepeat > 10) { + std::cout << "WARNING Average ray repeated " << avgRayRepeat <<" times!!" << std::endl << "Consider lower sample to seed ratio" << std::endl; } } @@ -134,7 +136,8 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu int mutation_counter = 0; if (mutation_number == 0) - inputParticle.weight = b / inputParticle.probability.totalProbability(); + if (seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(); //poprawne dla próbek po wyborze z puli + else inputParticle.weight = 1; //poprawne dla losowych probek while (mutation_counter < mutation_number) { @@ -146,18 +149,29 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); - if (rejectionSampling) { + //// The if block calculate sample weights for diferent cases + // Treat sample that is not large step + if (!mutatedParticle.largeStep && rejectionSampling){ // Rejection sampling inputParticle.weight += b / mutation_number / inputParticle.probability.totalProbability() * (1 - ratio); mutatedParticle.weight += b / mutation_number / mutatedParticle.probability.totalProbability() * ratio; - } - else { + } + else if(!mutatedParticle.largeStep){ // No rejection sampling inputParticle.weight = b / inputParticle.probability.totalProbability(); mutatedParticle.weight = b / mutatedParticle.probability.totalProbability(); + } + // Treat sample that is a large step (other cases treated before), with probability of generation is equal to 1 + else if (rejectionSampling){ + // Large step, rejection sampling + mutatedParticle.weight = 1. / mutation_number; + } + else { + // Large step, no rejection sampling + mutatedParticle.weight = 1; } - + //// Make particle seletion, update and SR casting if(mutatedParticle.largeStep==1) { if (rejectionSampling) CastShadowRays(inputParticle); @@ -245,7 +259,7 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML mutetedParticle.travelProbability.clear(); mutetedParticle.matAbsorbtions.clear(); mutetedParticle.largeStep = 1; - mutetedParticle.weight = 0; + mutetedParticle.weight =0; mutetedParticle.sourceDir1 = GetRandValue(); mutetedParticle.sourceDir2 = GetRandValue(); diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 038bcdab..c849b7b9 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -212,7 +212,6 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); - largeStepProb = 0; // to be deleted after proper large step implementation, currently fixed to avoid errors std::cout << "Mutation number: " << mutationNumber << std::endl; std::cout << "Large step probability: " << largeStepProb << std::endl << "Currently large steps are disabled/not implemented" << std::endl; From 2f875e1717f950e173cb22e63329c48063414840 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 18:40:01 +0200 Subject: [PATCH 13/56] Working rejection sampling --- src/spps-agh/calculation_cores/MLTCore.cpp | 30 ++++++---------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 4bd725f7..6718b914 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -151,34 +151,20 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu //// The if block calculate sample weights for diferent cases // Treat sample that is not large step - if (!mutatedParticle.largeStep && rejectionSampling){ + if (rejectionSampling){ // Rejection sampling - inputParticle.weight += b / mutation_number / inputParticle.probability.totalProbability() * (1 - ratio); - mutatedParticle.weight += b / mutation_number / mutatedParticle.probability.totalProbability() * ratio; + inputParticle.weight += (b / inputParticle.probability.totalProbability() * (1 - ratio)) / (mutation_number * (1 + pLarge)) ; + mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability() * ratio) + mutatedParticle.largeStep) / (mutation_number * (1 + pLarge)); } - else if(!mutatedParticle.largeStep){ + else{ // No rejection sampling - inputParticle.weight = b / inputParticle.probability.totalProbability(); - mutatedParticle.weight = b / mutatedParticle.probability.totalProbability(); - } - // Treat sample that is a large step (other cases treated before), with probability of generation is equal to 1 - else if (rejectionSampling){ - // Large step, rejection sampling - mutatedParticle.weight = 1. / mutation_number; - } - else { - // Large step, no rejection sampling - mutatedParticle.weight = 1; + inputParticle.weight = b / inputParticle.probability.totalProbability() / (1 + pLarge); + mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability() + mutatedParticle.largeStep) / (1 + pLarge); } //// Make particle seletion, update and SR casting - if(mutatedParticle.largeStep==1) - { - if (rejectionSampling) CastShadowRays(inputParticle); - inputParticle = mutatedParticle; - inputParticle.largeStep = 0; - } - else if(ratio > GetRandValue()) + mutatedParticle.largeStep = 0; + if(ratio > GetRandValue()) { if (rejectionSampling) CastShadowRays(inputParticle); inputParticle = mutatedParticle; From c8b81f11e42f13e78d0281f88a2fe36d1e89171c Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 20:38:18 +0200 Subject: [PATCH 14/56] Remove large step warning --- src/spps-agh/data_manager/core_configurationAGH.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index c849b7b9..dde554ba 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -213,8 +213,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); std::cout << "Mutation number: " << mutationNumber << std::endl; - std::cout << "Large step probability: " << largeStepProb << std::endl << - "Currently large steps are disabled/not implemented" << std::endl; + std::cout << "Large step probability: " << largeStepProb << std::endl; std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << "Currently spec. ref. multiplier is disabled" << std::endl; std::cout << "Refection penalty multiplier :" << reflPenaltyMultip << std::endl; From 1fa254f671ff4e295d89c7b7263b6aad84e02e07 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 20:39:17 +0200 Subject: [PATCH 15/56] Increment SPPS-AGH verion number to 0.4.0 --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index d97460d0..6b5be8f4 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -31,7 +31,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.36 (16.06.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.0 (16.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From b8cc8102f3f298ef1d46a4efff61d3d90b9fecdb Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 16 Oct 2021 22:56:48 +0200 Subject: [PATCH 16/56] Fix MLT for 0 mutations --- src/spps-agh/calculation_cores/MLTCore.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 6718b914..764267a9 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -136,7 +136,7 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu int mutation_counter = 0; if (mutation_number == 0) - if (seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(); //poprawne dla próbek po wyborze z puli + if (!seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(); //poprawne dla próbek po wyborze z puli else inputParticle.weight = 1; //poprawne dla losowych probek while (mutation_counter < mutation_number) From 2e88a17830a94bf8c3272629395259b0d99d701c Mon Sep 17 00:00:00 2001 From: wbinek Date: Sun, 17 Oct 2021 00:41:58 +0200 Subject: [PATCH 17/56] Increment SPPS-AGH version to 0.4.1 --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 6b5be8f4..5d0c1c36 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -31,7 +31,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.0 (16.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.1 (17.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 4053f311581a9273ffb669dd7f28d5bc2fbe44eb Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 19 Oct 2021 12:17:40 +0200 Subject: [PATCH 18/56] Add separate variable for enabling seed selection --- .../tree_core/e_core_spps_aghcore_advanced.h | 13 +++++++++++-- src/spps-agh/calculation_cores/MLTCore.cpp | 7 ++++--- src/spps-agh/calculation_cores/MLTCore.h | 1 + .../data_manager/core_configurationAGH.cpp | 14 ++++++++------ src/spps-agh/data_manager/core_configurationAGH.h | 3 ++- src/spps-agh/sppsNeeAGH.cpp | 3 ++- 6 files changed, 28 insertions(+), 13 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index c050494b..a83b26a8 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -251,6 +251,14 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element void Modified(Element* eModif) { + t_elementInfo filsInfo = eModif->GetElementInfos(); + if (filsInfo.libelleElement == "selectSeeds") { + this->SetReadOnlyConfig("seedSampleRatio", !this->GetBoolConfig("selectSeeds")); + } + else if (filsInfo.libelleElement == "mutationNumber") { + this->SetReadOnlyConfig("largeStepProb", !(bool)this->GetIntegerConfig("mutationNumber")); + } + Element::Modified(eModif); } @@ -275,10 +283,11 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element } void InitAdditionalMLTProperties(Element* confCore) { confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use - Rejection sampling"), true, true); - confCore->AppendPropertyBool("useFeatureScaling", wxTRANSLATE("Use - Feature scaling"), true, true); + confCore->AppendPropertyBool("useWeightDecayCompensation", wxTRANSLATE("Use - Weight decay compensation"), true, true); confCore->AppendPropertyBool("mutateSource", wxTRANSLATE("Use - Mutate source"), true, true); + confCore->AppendPropertyBool("selectSeeds", wxTRANSLATE("Use - Seed selection"), true, true); - confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0.01, true); confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 764267a9..88a61dd9 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -17,6 +17,7 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT reflectionPenaltyMultip = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_REFL_PENALTY_MULTIP); seedSampleRatio = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); + selectSeeds = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_SELECT_SEEDS); rejectionSampling = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_REJECTION_SAMPLING); costMatAbsorption = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_ABSORPTION); costMatDiffusion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_MAT_DIFFUSION); @@ -57,9 +58,9 @@ void MLTCore::ChooseSeeds() { InitialSeedsSize = InitialSeeds.size(); // If seed chosing disabled - if (seedSampleRatio == 0) { + if (!selectSeeds) { Seeds = InitialSeeds; - std::cout << "Seed choosing disabled, copied all initial seeds" << std::endl; + std::cout << "Seed selection disabled, copied all initial seeds" << std::endl; return; } @@ -123,7 +124,7 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) input_particle.shadowRays.pop_front(); //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; - if (!seedSampleRatio == 0) shadowRay.energie *= InitialSeedsSize / SeedsSize; // Compensate energy for seed selection + if (selectSeeds) shadowRay.energie *= InitialSeedsSize / SeedsSize; // Compensate energy for seed selection Run(shadowRay); //applicationTools.outputTool->SaveParticule(); } diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 80350e19..d7315b27 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -21,6 +21,7 @@ class MLTCore : public NextEventEstimationCore float reflectionPenaltyMultip = 0.95; float seedSampleRatio = 0; + bool selectSeeds = true; bool rejectionSampling = true; bool costMatAbsorption = true; bool costMatDiffusion = false; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index dde554ba..3a27733d 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -209,6 +209,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int costShadowRays = advancedNode->GetProperty("costShadowRays").ToInt(); int costAirAbsoprtion = advancedNode->GetProperty("costAirAbs").ToInt(); int mutateSource = advancedNode->GetProperty("mutateSource").ToInt(); + int selectSeeds = advancedNode->GetProperty("selectSeeds").ToInt(); float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); @@ -217,9 +218,9 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << "Currently spec. ref. multiplier is disabled" << std::endl; std::cout << "Refection penalty multiplier :" << reflPenaltyMultip << std::endl; - std::cout << "Seed to sample rato :" << seedSampleRatio << std::endl; - if(seedSampleRatio==0) std::cout << "Seed selection disabled!" << std::endl; - + std::cout << "Seed selection: " << selectSeeds << std::endl; + if (selectSeeds) std::cout << "Seed to sample rato :" << seedSampleRatio << std::endl; + std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; std::cout << "Rejecton sampling: " << rejectionSampling << std::endl; std::cout << "Cost - absorption: " << costMatAbsorption << std::endl; @@ -228,8 +229,8 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Cost - reflection dir.: " << costReflDir << std::endl; std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; - if (advancedNode->GetProperty("useFeatureScaling").ToInt()) { - std::cout << "Feature scaling not implemented yet!!" << std::endl; + if (advancedNode->GetProperty("useWeightDecayCompensation").ToInt()) { + std::cout << "Weight Decay Compensation not implemented yet!!" << std::endl; } SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); @@ -238,8 +239,9 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(FPROP_MLT_REFL_PENALTY_MULTIP, reflPenaltyMultip); SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, rejectionSampling); - SetConfigInformation(IPROP_MLT_DO_FEATURE_SCALING, advancedNode->GetProperty("useFeatureScaling").ToInt()); + SetConfigInformation(IPROP_MLT_DO_WEIGHT_DECAY_COMPENSATION, advancedNode->GetProperty("useWeightDecayCompensation").ToInt()); SetConfigInformation(IPROP_MLT_DO_MUTATE_SOURCE, mutateSource); + SetConfigInformation(IPROP_MLT_DO_SELECT_SEEDS, selectSeeds); SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, costMatAbsorption); SetConfigInformation(IPROP_MLT_COST_MAT_DIFFUSION, costMatDiffusion); SetConfigInformation(IPROP_MLT_COST_SHADOW_RAYS, costShadowRays); diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index b6509dd4..9d647d92 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -45,8 +45,9 @@ class Core_ConfigurationAGH : public Core_Configuration //MLT related boolean properties IPROP_MLT_DO_REJECTION_SAMPLING, - IPROP_MLT_DO_FEATURE_SCALING, + IPROP_MLT_DO_WEIGHT_DECAY_COMPENSATION, IPROP_MLT_DO_MUTATE_SOURCE, + IPROP_MLT_DO_SELECT_SEEDS, IPROP_MLT_COST_MAT_ABSORPTION, IPROP_MLT_COST_MAT_DIFFUSION, IPROP_MLT_COST_SHADOW_RAYS, diff --git a/src/spps-agh/sppsNeeAGH.cpp b/src/spps-agh/sppsNeeAGH.cpp index ab2d12fc..fbb73994 100644 --- a/src/spps-agh/sppsNeeAGH.cpp +++ b/src/spps-agh/sppsNeeAGH.cpp @@ -91,7 +91,8 @@ void runSourceCalculation( progressOperation* parentOperation, t_ToolBox& applic int progresSteps = quandparticules; if (*applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::I_PROP_CALCULATION_CORE_SELLECTION) == MLT) { int mutationNumber = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); - float seedSampleRatio = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); + bool selectSeeds = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_SELECT_SEEDS); + float seedSampleRatio = (selectSeeds) ? *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO) : 1; progresSteps += seedSampleRatio*quandparticules* mutationNumber; } From b282fd7efe21f5587068d0831daf75f781d89bfc Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 19 Oct 2021 12:29:08 +0200 Subject: [PATCH 19/56] Add mutation generator parameters as input --- .../tree_core/e_core_spps_aghcore_advanced.h | 2 ++ src/spps-agh/calculation_cores/MLTCore.cpp | 5 ++++- src/spps-agh/calculation_cores/MLTCore.h | 3 +++ src/spps-agh/data_manager/core_configurationAGH.cpp | 10 +++++++++- src/spps-agh/data_manager/core_configurationAGH.h | 4 +++- 5 files changed, 21 insertions(+), 3 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index a83b26a8..13e02233 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -289,6 +289,8 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0.01, true); confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); + confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); + confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 88a61dd9..92fd2168 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -25,6 +25,9 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT costShadowRays = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_SHADOW_RAYS); costAirAbsoprtion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_AIR_ABSORPTION); mutateSource = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_MUTATE_SOURCE); + + genS1 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S1); + genS2 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S2); }; bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) @@ -497,7 +500,7 @@ bool MLTCore::GetNextCollision(CONF_PARTICULE_MLT& configurationP, double& dista void MLTCore::MutationGenerator(double& value) const { - float s1 = 1. / 1024, s2 = 1. / 64; + float s1 = 1. / genS1, s2 = 1. / genS2; float dv = s2*exp(-log(s2 / s1)*GetRandValue()); if (GetRandValue() < 0.5) { value += dv; if (value > 1) value -= 1; diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index d7315b27..046be213 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -30,6 +30,9 @@ class MLTCore : public NextEventEstimationCore bool costAirAbsoprtion = false; bool mutateSource = true; + int genS1 = 1024; + int genS2 = 64; + protected: std::list InitialSeeds; std::list Seeds; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 3a27733d..2fb984cc 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -213,6 +213,9 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); + int mutationGeneratorS1 = advancedNode->GetProperty("mutationGeneratorS1").ToInt(); + int mutationGeneratorS2 = advancedNode->GetProperty("mutationGeneratorS2").ToInt(); + std::cout << "Mutation number: " << mutationNumber << std::endl; std::cout << "Large step probability: " << largeStepProb << std::endl; std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << @@ -228,10 +231,12 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; std::cout << "Cost - reflection dir.: " << costReflDir << std::endl; std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; - if (advancedNode->GetProperty("useWeightDecayCompensation").ToInt()) { std::cout << "Weight Decay Compensation not implemented yet!!" << std::endl; } + + std::cout << "Mutation generator S1: " << mutationGeneratorS1 << std::endl; + std::cout << "Mutation generator S2: " << mutationGeneratorS2 << std::endl; SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, largeStepProb); @@ -248,6 +253,9 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_COST_REFL_DIR, costReflDir); SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, costAirAbsoprtion); SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, seedSampleRatio); + + SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S1, mutationGeneratorS1); + SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S2, mutationGeneratorS2); } } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index 9d647d92..5d5db2be 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -52,7 +52,9 @@ class Core_ConfigurationAGH : public Core_Configuration IPROP_MLT_COST_MAT_DIFFUSION, IPROP_MLT_COST_SHADOW_RAYS, IPROP_MLT_COST_REFL_DIR, - IPROP_MLT_COST_AIR_ABSORPTION + IPROP_MLT_COST_AIR_ABSORPTION, + IPROP_MLT_MUTATION_GENERATOR_S1, + IPROP_MLT_MUTATION_GENERATOR_S2, }; /** * Initialisation des paramètres du coeur de calcul à partir d'un fichier XML From e51dff783da0a12c08c49c9946e7d4b88f94e9df Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 19 Oct 2021 19:47:26 +0200 Subject: [PATCH 20/56] Add initial take on weight decay compensation --- .../tree_core/e_core_spps_aghcore_advanced.h | 2 +- src/spps-agh/calculation_cores/MLTCore.cpp | 59 +++++++++++++------ src/spps-agh/calculation_cores/MLTCore.h | 6 +- .../data_manager/core_configurationAGH.cpp | 12 ++-- src/spps-agh/sppsNeeAGH.cpp | 2 + src/spps-agh/sppsNeeAGHTypes.h | 11 +++- 6 files changed, 64 insertions(+), 28 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 13e02233..3563ca39 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -288,7 +288,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyBool("selectSeeds", wxTRANSLATE("Use - Seed selection"), true, true); confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0.01, true); - confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 2, 0.01, true); + confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 92fd2168..a53bee81 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -25,6 +25,7 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT costShadowRays = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_SHADOW_RAYS); costAirAbsoprtion = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_COST_AIR_ABSORPTION); mutateSource = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_MUTATE_SOURCE); + weightDecayCompensation = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_WEIGHT_DECAY_COMPENSATION); genS1 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S1); genS2 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S2); @@ -47,7 +48,8 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) RunInitialSeed(inputParticle); InitialSeeds.push_back(inputParticle); - b += inputParticle.probability.totalProbability() / totalSeedNumber; + + //b += inputParticle.probability.totalProbability() / totalSeedNumber; } else { @@ -57,6 +59,25 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) return true; } +void MLTCore::CalculatewDC() { + if (weightDecayCompensation) + { + int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); + wDC = 0; + for (CONF_PARTICULE_MLT seed : InitialSeeds) { + double currentProbabilityRaw = seed.probability.currentProbability / pow(reflectionPenaltyMultip, seed.reflectionOrder); + wDC += pow(currentProbabilityRaw, (1. / (seed.reflectionsNum+1))) / totalSeedNumber; + } + } +} + +void MLTCore::CalculateB() { + int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); + for (CONF_PARTICULE_MLT seed : InitialSeeds) { + b += seed.probability.totalProbability(wDC) / totalSeedNumber; + } +} + void MLTCore::ChooseSeeds() { InitialSeedsSize = InitialSeeds.size(); @@ -76,7 +97,7 @@ void MLTCore::ChooseSeeds() { for (CONF_PARTICULE_MLT seed : InitialSeeds) { accounted = false; - float currentScore = seed.probability.totalProbability() / b * (N_target/totalSeedNumber); + float currentScore = seed.probability.totalProbability(wDC) / b * (N_target/totalSeedNumber); if (currentScore > 1) { individualRays++; accounted = true; @@ -114,6 +135,7 @@ void MLTCore::RunMutation() { CONF_PARTICULE_MLT inputParticle = Seeds.front(); Seeds.pop_front(); + inputParticle.isInitialSeed = false; DoMutations(inputParticle, (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size())); } @@ -140,7 +162,7 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu int mutation_counter = 0; if (mutation_number == 0) - if (!seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(); //poprawne dla próbek po wyborze z puli + if (!seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(wDC); //poprawne dla próbek po wyborze z puli else inputParticle.weight = 1; //poprawne dla losowych probek while (mutation_counter < mutation_number) @@ -148,7 +170,7 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu Mutate(inputParticle, mutatedParticle); Follow(mutatedParticle); - ratio = std::min(mutatedParticle.probability.totalProbability() / inputParticle.probability.totalProbability(), 1.); + ratio = std::min(mutatedParticle.probability.totalProbability(wDC) / inputParticle.probability.totalProbability(wDC), 1.); //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); @@ -157,13 +179,13 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu // Treat sample that is not large step if (rejectionSampling){ // Rejection sampling - inputParticle.weight += (b / inputParticle.probability.totalProbability() * (1 - ratio)) / (mutation_number * (1 + pLarge)) ; - mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability() * ratio) + mutatedParticle.largeStep) / (mutation_number * (1 + pLarge)); + inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1 - ratio)) / (mutation_number * (1 + pLarge)) ; + mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio) + mutatedParticle.largeStep) / (mutation_number * (1 + pLarge)); } else{ // No rejection sampling - inputParticle.weight = b / inputParticle.probability.totalProbability() / (1 + pLarge); - mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability() + mutatedParticle.largeStep) / (1 + pLarge); + inputParticle.weight = b / inputParticle.probability.totalProbability(wDC) / (1 + pLarge); + mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability(wDC) + mutatedParticle.largeStep) / (1 + pLarge); } //// Make particle seletion, update and SR casting @@ -190,6 +212,7 @@ bool MLTCore::RunInitialSeed(CONF_PARTICULE_MLT& inputParticle) ParticleDistribution::GenSphereDistribution(inputParticle, distanceToTravel, &rndDir1, &rndDir2); inputParticle.sourceDir1=rndDir1; inputParticle.sourceDir2=rndDir2; + inputParticle.isInitialSeed=true; SetNextParticleCollision(inputParticle); //1er test de collision SetNextParticleCollisionWithObstructionElement(inputParticle); // Test de collision avec objet virtuel encombrant @@ -238,6 +261,7 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML mutetedParticle = input_particle; mutetedParticle.shadowRays.clear(); mutetedParticle.probability.partialProbability.clear(); + mutetedParticle.probability.srRelfectionOrder.clear(); mutetedParticle.weight = 0; //bool large_step = (GetRandValue() < pLarge) ? 1 : 0; @@ -339,18 +363,17 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn std::vector srProbabilities; GenerateShadowRays(configurationP, materialInfo, faceNormal, deltaT, distanceToTravel, configurationP.shadowRays, faceInfo, &srProbabilities); - if (configurationP.reflectionOrder > 0) + if(costShadowRays) // If shadow ray cost is enabled adds the currentProb * srProb for each Shadow Ray { - if(costShadowRays) // If shadow ray cost is enabled adds the currentProb * srProb for each Shadow Ray - { - for (double srProb : srProbabilities) { - configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability * srProb); - } - } - else { // If shadow ray cost is disablerd adds the currentProb once (this option is rather useless) - configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability); + for (double srProb : srProbabilities) { + configurationP.probability.srRelfectionOrder.push_back(configurationP.reflectionOrder); + configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability * srProb); } } + else { // If shadow ray cost is disabled adds the currentProb once (this option is rather useless) + configurationP.probability.srRelfectionOrder.push_back(configurationP.reflectionOrder); + configurationP.probability.partialProbability.push_back(configurationP.probability.currentProbability); + } //******************************* REFLECTION ***************************************// @@ -380,7 +403,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn if(costReflDir) configurationP.probability.currentProbability *= dirProbability; configurationP.probability.currentProbability *= reflectionPenaltyMultip; - + return false; } diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 046be213..9991b909 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -12,8 +12,11 @@ class MLTCore : public NextEventEstimationCore bool SeedIsEmpty(); void RunMutation(); void ChooseSeeds(); + void CalculateB(); + void CalculatewDC(); - float b = 0; //bias value - auto calculated - start with 0!; + float b = 0; //bias value - auto calculated - start with 0!; + float wDC = 1; //weight decay compensation factor (auto set if option enabled); float pLarge = 0.5; float pSpecular = 0.5; @@ -29,6 +32,7 @@ class MLTCore : public NextEventEstimationCore bool costShadowRays = true; bool costAirAbsoprtion = false; bool mutateSource = true; + bool weightDecayCompensation = true; int genS1 = 1024; int genS2 = 64; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 2fb984cc..7d8dc65f 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -212,6 +212,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int selectSeeds = advancedNode->GetProperty("selectSeeds").ToInt(); float seedSampleRatio = advancedNode->GetProperty("seedSampleRatio").ToFloat(); int rejectionSampling = advancedNode->GetProperty("useRejectionSampling").ToInt(); + int useWeightDecayCompensation = advancedNode->GetProperty("useWeightDecayCompensation").ToInt(); int mutationGeneratorS1 = advancedNode->GetProperty("mutationGeneratorS1").ToInt(); int mutationGeneratorS2 = advancedNode->GetProperty("mutationGeneratorS2").ToInt(); @@ -224,6 +225,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Seed selection: " << selectSeeds << std::endl; if (selectSeeds) std::cout << "Seed to sample rato :" << seedSampleRatio << std::endl; + std::cout << "Weight Decay Compensation: " << useWeightDecayCompensation << std::endl; std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; std::cout << "Rejecton sampling: " << rejectionSampling << std::endl; std::cout << "Cost - absorption: " << costMatAbsorption << std::endl; @@ -231,12 +233,12 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Cost - ShadowRays: " << costShadowRays << std::endl; std::cout << "Cost - reflection dir.: " << costReflDir << std::endl; std::cout << "Cost - air absoprtion: " << costAirAbsoprtion << std::endl; - if (advancedNode->GetProperty("useWeightDecayCompensation").ToInt()) { - std::cout << "Weight Decay Compensation not implemented yet!!" << std::endl; - } - + std::cout << "Mutation generator S1: " << mutationGeneratorS1 << std::endl; std::cout << "Mutation generator S2: " << mutationGeneratorS2 << std::endl; + + if (useWeightDecayCompensation && reflPenaltyMultip > 0.95) + std::cout << "WARNING! Using Weight Decay Compensation with high relfection penalty multiplier may cause artefacts!" << std::endl; SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, largeStepProb); @@ -244,7 +246,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(FPROP_MLT_REFL_PENALTY_MULTIP, reflPenaltyMultip); SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, rejectionSampling); - SetConfigInformation(IPROP_MLT_DO_WEIGHT_DECAY_COMPENSATION, advancedNode->GetProperty("useWeightDecayCompensation").ToInt()); + SetConfigInformation(IPROP_MLT_DO_WEIGHT_DECAY_COMPENSATION, useWeightDecayCompensation); SetConfigInformation(IPROP_MLT_DO_MUTATE_SOURCE, mutateSource); SetConfigInformation(IPROP_MLT_DO_SELECT_SEEDS, selectSeeds); SetConfigInformation(IPROP_MLT_COST_MAT_ABSORPTION, costMatAbsorption); diff --git a/src/spps-agh/sppsNeeAGH.cpp b/src/spps-agh/sppsNeeAGH.cpp index fbb73994..2e9637cd 100644 --- a/src/spps-agh/sppsNeeAGH.cpp +++ b/src/spps-agh/sppsNeeAGH.cpp @@ -177,6 +177,8 @@ void runSourceCalculation( progressOperation* parentOperation, t_ToolBox& applic if(*applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::I_PROP_CALCULATION_CORE_SELLECTION)==MLT) { int mutationNumber = *applicationTools.configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); + static_cast(applicationTools.calculationTool)->CalculatewDC(); + static_cast(applicationTools.calculationTool)->CalculateB(); static_cast(applicationTools.calculationTool)->ChooseSeeds(); // choose seeds from the initial seed pool while (!static_cast(applicationTools.calculationTool)->SeedIsEmpty()) { diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 5d0c1c36..71345472 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -31,7 +31,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.1 (17.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.2 (19.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION @@ -69,17 +69,22 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH int reflectionsNum = 0; int creationTime = 0; int largeStep = 0; + bool isInitialSeed = false; class { public: double currentProbability = 1; //probability of generating ray path up to some point + double currentProbabilityRaw = 1; //probability of generating ray payth up to some point without correction std::vector partialProbability; //probability of generating shadow ray + std::vector srRelfectionOrder; //Used for weight decay compensation - double totalProbability() //sum of shadow rays probability - total ray contribution + double totalProbability(double wDC) //sum of shadow rays probability - total ray contribution { double totalProb = 0; - for(double val : partialProbability){ + for (int i = 0; i < partialProbability.size();i++) { + double val = partialProbability[i]; + val = val / pow(wDC, srRelfectionOrder[i]); totalProb += val; } return totalProb; From ac5728a6b7ffab003d75639838b07922ea87d596 Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 19 Oct 2021 23:52:58 +0200 Subject: [PATCH 21/56] Change default parameter value --- .../data_manager/tree_core/e_core_spps_aghcore_advanced.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 3563ca39..610b8bcf 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -288,7 +288,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyBool("selectSeeds", wxTRANSLATE("Use - Seed selection"), true, true); confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0.01, true); - confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.95, false, 2, true, true, 10, 0.01, true); + confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.8, false, 2, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); From 8674edd49248cca0f61f61029e84ad9dde51f91c Mon Sep 17 00:00:00 2001 From: wbinek Date: Wed, 20 Oct 2021 00:00:15 +0200 Subject: [PATCH 22/56] Fix regression, seed selection without mutations error --- src/spps-agh/calculation_cores/MLTCore.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index a53bee81..2e01578f 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -162,7 +162,7 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu int mutation_counter = 0; if (mutation_number == 0) - if (!seedSampleRatio==0) inputParticle.weight = b / inputParticle.probability.totalProbability(wDC); //poprawne dla próbek po wyborze z puli + if (selectSeeds) inputParticle.weight = b / inputParticle.probability.totalProbability(wDC); //poprawne dla próbek po wyborze z puli else inputParticle.weight = 1; //poprawne dla losowych probek while (mutation_counter < mutation_number) From 97336777a0578a8b7a9637f02467c7aaddaaeb86 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 22 Oct 2021 13:15:40 +0200 Subject: [PATCH 23/56] Fix very rare NaN in wDC calculation --- src/spps-agh/calculation_cores/MLTCore.cpp | 5 ++++- src/spps-agh/tools/dotreflectionAGH.h | 4 +--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 2e01578f..777fa320 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -66,8 +66,10 @@ void MLTCore::CalculatewDC() { wDC = 0; for (CONF_PARTICULE_MLT seed : InitialSeeds) { double currentProbabilityRaw = seed.probability.currentProbability / pow(reflectionPenaltyMultip, seed.reflectionOrder); - wDC += pow(currentProbabilityRaw, (1. / (seed.reflectionsNum+1))) / totalSeedNumber; + double wDCIncrement = pow(currentProbabilityRaw, (1. / (seed.reflectionsNum + 1))) / totalSeedNumber; + wDC += wDCIncrement; } + if (isnan(wDC)) std::cout << "!!!!!! wDC is NaN !!!!!!!!!" << std::endl; } } @@ -76,6 +78,7 @@ void MLTCore::CalculateB() { for (CONF_PARTICULE_MLT seed : InitialSeeds) { b += seed.probability.totalProbability(wDC) / totalSeedNumber; } + if (isnan(b)) std::cout << "!!!!!! b is NaN !!!!!!!!!" << std::endl; } void MLTCore::ChooseSeeds() { diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index db818579..be4e21db 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -66,8 +66,6 @@ class ReflectionLawsAGH : public ReflectionLaws target = rotationMatrix*target; target.normalize(); - double test = target.dot(faceNormal); - //test if reflected ray is pointing into ground if(target.dot(faceNormal) <= 0) { @@ -168,7 +166,7 @@ class ReflectionLawsMLT : public ReflectionLawsAGH //probability = ((expo + 2) / (2 * M_PI)) * pow(cos(phi),expo); //probability = (1/M_2PI) * pow(cos(phi), expo); probability = pow(cos(phi), expo); //normalization factor excluded, to aviod large numbers - if (isnan(probability)) [[unlikely]] { + if (isnan(probability) || probability<0) [[unlikely]] { probability = 0; } From 1e1957ca2b3e316ee2b8e5c5d72bb35a9e685d80 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 22 Oct 2021 16:39:55 +0200 Subject: [PATCH 24/56] Increment SPPS-AGH version number --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 71345472..51e96ff8 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -31,7 +31,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.2 (19.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.3 (22.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 147f93a6ef788a62698f9538864cabe8831c8cb5 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 22 Oct 2021 22:51:15 +0200 Subject: [PATCH 25/56] Potential fix for Phong normalization regression --- src/spps-agh/tools/brdfreflection.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/spps-agh/tools/brdfreflection.cpp b/src/spps-agh/tools/brdfreflection.cpp index b24a3538..2668c14a 100644 --- a/src/spps-agh/tools/brdfreflection.cpp +++ b/src/spps-agh/tools/brdfreflection.cpp @@ -42,7 +42,7 @@ int RaySphereIntersection(const vec3& p1, const vec3& p2, const vec3& sc, double double BRDFs::calculatePhongNormalizationFactor(const vec3& faceNormal, const vec3& specularDirection, const int& n) { - double d = std::min((float)faceNormal.dot(specularDirection),1.f); + double d = std::max(std::min((float)faceNormal.dot(specularDirection),1.f), 0.f); double PNF; //double PNF = BRDFs::calculatePNF(d, n); if (n > 1000) [[unlikely]] { @@ -50,7 +50,12 @@ double BRDFs::calculatePhongNormalizationFactor(const vec3& faceNormal, const ve } else { int d_idx = std::round(d * 1000); - PNF = BRDFs::phongFactorTable[n][d_idx]; + try { + PNF = BRDFs::phongFactorTable[n][d_idx]; \ + } + catch(const std::exception& e){ + std::cout << e.what() << std::endl; + } } return PNF; From b03f0fce8f903a13d1ebd1d46d64ba3b29598f34 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 23 Oct 2021 21:24:47 +0200 Subject: [PATCH 26/56] Release memory when seeds are selected --- src/spps-agh/calculation_cores/MLTCore.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 777fa320..1f9f970d 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -127,6 +127,8 @@ void MLTCore::ChooseSeeds() { std::cout << "WARNING Average ray repeated " << avgRayRepeat <<" times!!" << std::endl << "Consider lower sample to seed ratio" << std::endl; } + // Release memory taken by initial seeds, as they are no more reqired + InitialSeeds.clear(); } bool MLTCore::SeedIsEmpty() From a542f17989dafb8e620bb2deda9c5812b5d02379 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 23 Oct 2021 21:25:17 +0200 Subject: [PATCH 27/56] Remove try catch clause used for debug purposes --- src/spps-agh/tools/brdfreflection.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/spps-agh/tools/brdfreflection.cpp b/src/spps-agh/tools/brdfreflection.cpp index 2668c14a..ea11486e 100644 --- a/src/spps-agh/tools/brdfreflection.cpp +++ b/src/spps-agh/tools/brdfreflection.cpp @@ -50,12 +50,7 @@ double BRDFs::calculatePhongNormalizationFactor(const vec3& faceNormal, const ve } else { int d_idx = std::round(d * 1000); - try { - PNF = BRDFs::phongFactorTable[n][d_idx]; \ - } - catch(const std::exception& e){ - std::cout << e.what() << std::endl; - } + PNF = BRDFs::phongFactorTable[n][d_idx]; } return PNF; From 64c29697002a891375eb95e4336a9a5774ce0b13 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 23 Oct 2021 21:25:52 +0200 Subject: [PATCH 28/56] Add propability for default case in dotreflection --- src/spps-agh/tools/dotreflectionAGH.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index be4e21db..0b2381da 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -114,6 +114,7 @@ class ReflectionLawsMLT : public ReflectionLawsAGH case REFLECTION_LAW_PHONG: return BaseWnReflection(vectorDirection, faceNormal, 1, rnd1, rnd2, probability); default: + probability = 1; return SpecularReflection(vectorDirection, faceNormal); }; } From 7eea9f4c42c6107447498efa9f16ef30029dc801 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 23 Oct 2021 21:54:24 +0200 Subject: [PATCH 29/56] Add some additional debug comunicates --- src/spps-agh/calculation_cores/MLTCore.cpp | 5 +---- src/spps-agh/sppsNeeAGHTypes.h | 1 + 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 1f9f970d..fda954c3 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -44,12 +44,8 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) inputParticle.SetInitialState(inputParticle); inputParticle.weight = 0; - int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); - RunInitialSeed(inputParticle); InitialSeeds.push_back(inputParticle); - - //b += inputParticle.probability.totalProbability() / totalSeedNumber; } else { @@ -177,6 +173,7 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu ratio = std::min(mutatedParticle.probability.totalProbability(wDC) / inputParticle.probability.totalProbability(wDC), 1.); + if (isnan(ratio)) [[unlikely]] std::cout << "Ratio is NaN, may cause serious erros" << std::endl; //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 51e96ff8..61bcdf0a 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -87,6 +87,7 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH val = val / pow(wDC, srRelfectionOrder[i]); totalProb += val; } + if (isnan(totalProb)) [[unlikely]] std::cout << "Total probability of particle is NaN, this may lead to errors" << std::endl; return totalProb; } From ac27b610f35bcbdeaa7d7392cc8066ef1d119c11 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sun, 24 Oct 2021 11:33:56 +0200 Subject: [PATCH 30/56] Include iostream to enable cout for PARTICULE_MLT --- src/spps-agh/sppsNeeAGHTypes.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 61bcdf0a..04490b94 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -25,6 +25,7 @@ //#include "coreTypes.h" #include "sppsTypes.h" #include +#include #ifndef SPPSNEE_AGH_TYPES #define SPPSNEE_AGH_TYPES From 9447893636b678305a28277142f22ad1c504f407 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sun, 24 Oct 2021 11:37:02 +0200 Subject: [PATCH 31/56] Add new parameter - seed selection exponent, disable specular_refl_prob parameter --- .../tree_core/e_core_spps_aghcore_advanced.h | 4 +++- src/spps-agh/calculation_cores/MLTCore.cpp | 14 ++++++++++---- src/spps-agh/calculation_cores/MLTCore.h | 3 ++- .../data_manager/core_configurationAGH.cpp | 13 ++++++++----- src/spps-agh/data_manager/core_configurationAGH.h | 1 + 5 files changed, 24 insertions(+), 11 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 610b8bcf..b2b68909 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -254,6 +254,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element t_elementInfo filsInfo = eModif->GetElementInfos(); if (filsInfo.libelleElement == "selectSeeds") { this->SetReadOnlyConfig("seedSampleRatio", !this->GetBoolConfig("selectSeeds")); + this->SetReadOnlyConfig("seedSelectionExponent", !this->GetBoolConfig("selectSeeds")); } else if (filsInfo.libelleElement == "mutationNumber") { this->SetReadOnlyConfig("largeStepProb", !(bool)this->GetIntegerConfig("mutationNumber")); @@ -278,7 +279,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element void InitBaseMLTProperties(Element* confCore) { confCore->AppendPropertyDecimal("largeStepProb", wxTRANSLATE("Tune - Large step probability"), 0.5, false, 2, true, true, 1, 0, true); - confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Tune - Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); + //confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Tune - Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Sampling - Mutation number"), 25, true, false, true, 0, 0); } void InitAdditionalMLTProperties(Element* confCore) { @@ -291,6 +292,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.8, false, 2, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); + confCore->AppendPropertyDecimal("seedSelectionExponent", wxTRANSLATE("Tune - Seed selection exponent"), 0.5, false, 2, true, true, 1, 0.01, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index fda954c3..c67dba81 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -12,10 +12,11 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT doDirectSoundCalculation = true; pLarge = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_LARGE_STEP_PROB); - pSpecular = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SPECULAR_REFL_PROB); + //pSpecular = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SPECULAR_REFL_PROB); mutation_number = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_NUMBER); reflectionPenaltyMultip = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_REFL_PENALTY_MULTIP); seedSampleRatio = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); + seedSelectionExponent = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SELECTION_EXPONENT); selectSeeds = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_SELECT_SEEDS); rejectionSampling = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_REJECTION_SAMPLING); @@ -89,14 +90,19 @@ void MLTCore::ChooseSeeds() { // Seed choosing algorithm int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); - float N_target = seedSampleRatio * totalSeedNumber; + double N_target = seedSampleRatio * totalSeedNumber; + double total_score = 0; int individualRays = 0; bool accounted = false; + for (CONF_PARTICULE_MLT seed : InitialSeeds) { + total_score += pow(seed.probability.totalProbability(wDC), seedSelectionExponent); + } + for (CONF_PARTICULE_MLT seed : InitialSeeds) { accounted = false; - float currentScore = seed.probability.totalProbability(wDC) / b * (N_target/totalSeedNumber); + float currentScore = pow(seed.probability.totalProbability(wDC), seedSelectionExponent) / total_score * (N_target); if (currentScore > 1) { individualRays++; accounted = true; @@ -387,7 +393,7 @@ bool MLTCore::MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rn { nouvDirection = ReflectionLawsMLT::SolveDiffusePart(configurationP.direction, *materialInfo, faceNormal, configurationP, rndDir1, rndDir2, dirProbability); if (costMatDiffusion) configurationP.probability.currentProbability *= materialInfo->diffusion; - //configurationP.energie /= (1-pSpecular) + //configurationP.energie *= pSpecular //dirProbability /= (1 - pSpecular); } else diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 9991b909..b10b6229 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -19,10 +19,11 @@ class MLTCore : public NextEventEstimationCore float wDC = 1; //weight decay compensation factor (auto set if option enabled); float pLarge = 0.5; - float pSpecular = 0.5; + float pSpecular = 1; int mutation_number = 25; float reflectionPenaltyMultip = 0.95; float seedSampleRatio = 0; + float seedSelectionExponent = 0.5; bool selectSeeds = true; bool rejectionSampling = true; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 7d8dc65f..bfd6203a 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -200,7 +200,8 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { //Load MLT parameters, do check on them, write to console if not implemented int mutationNumber = advancedNode->GetProperty("mutationNumber").ToInt(); float largeStepProb = advancedNode->GetProperty("largeStepProb").ToFloat(); - float specularReflProbabilityMulti = advancedNode->GetProperty("specularReflProbabilityMulti").ToFloat(); + //float specularReflProbabilityMulti = advancedNode->GetProperty("specularReflProbabilityMulti").ToFloat(); + float seedSelectionExponent = advancedNode->GetProperty("seedSelectionExponent").ToFloat(); float reflPenaltyMultip = advancedNode->GetProperty("reflectionPenalty").ToFloat(); int costMatAbsorption = advancedNode->GetProperty("costMatAbsorption").ToInt(); @@ -219,11 +220,12 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Mutation number: " << mutationNumber << std::endl; std::cout << "Large step probability: " << largeStepProb << std::endl; - std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << - "Currently spec. ref. multiplier is disabled" << std::endl; + //std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << + // "Currently spec. ref. multiplier is disabled" << std::endl; std::cout << "Refection penalty multiplier :" << reflPenaltyMultip << std::endl; std::cout << "Seed selection: " << selectSeeds << std::endl; - if (selectSeeds) std::cout << "Seed to sample rato :" << seedSampleRatio << std::endl; + if (selectSeeds) std::cout << " Seed to sample rato:" << seedSampleRatio << std::endl; + if (selectSeeds) std::cout << " Seed selection exponent:" << seedSelectionExponent << std::endl; std::cout << "Weight Decay Compensation: " << useWeightDecayCompensation << std::endl; std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; @@ -242,7 +244,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(FPROP_MLT_MUTATION_NUMBER, mutationNumber); SetConfigInformation(FPROP_MLT_LARGE_STEP_PROB, largeStepProb); - SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, specularReflProbabilityMulti); + //SetConfigInformation(FPROP_MLT_SPECULAR_REFL_PROB, specularReflProbabilityMulti); SetConfigInformation(FPROP_MLT_REFL_PENALTY_MULTIP, reflPenaltyMultip); SetConfigInformation(IPROP_MLT_DO_REJECTION_SAMPLING, rejectionSampling); @@ -255,6 +257,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_COST_REFL_DIR, costReflDir); SetConfigInformation(IPROP_MLT_COST_AIR_ABSORPTION, costAirAbsoprtion); SetConfigInformation(FPROP_MLT_SEED_SAMPLE_RATIO, seedSampleRatio); + SetConfigInformation(FPROP_MLT_SEED_SELECTION_EXPONENT, seedSelectionExponent); SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S1, mutationGeneratorS1); SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S2, mutationGeneratorS2); diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index 5d5db2be..39c05c78 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -32,6 +32,7 @@ class Core_ConfigurationAGH : public Core_Configuration FPROP_NEE_SHADOWRAY_PROB, FPROP_MLT_SEED_SAMPLE_RATIO, FPROP_MLT_REFL_PENALTY_MULTIP, + FPROP_MLT_SEED_SELECTION_EXPONENT, }; enum NEW_IPROP From dd5f5988f76253dc92ca8ac34a159edb348e467b Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 26 Oct 2021 12:47:57 +0200 Subject: [PATCH 32/56] Fix occasional crash on visability test --- src/spps-agh/CalculationCoreAGH.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spps-agh/CalculationCoreAGH.cpp b/src/spps-agh/CalculationCoreAGH.cpp index a4bc5f53..f218be76 100644 --- a/src/spps-agh/CalculationCoreAGH.cpp +++ b/src/spps-agh/CalculationCoreAGH.cpp @@ -225,6 +225,7 @@ bool CalculationCoreSPPS::VisabilityTest(const CONF_PARTICULE_AGH &configuration { testParticle.currentTetra = nextTetra; testParticle.nextModelIntersection.idface = GetTetraFaceCollision(testParticle, testParticle.direction, t); + if (testParticle.nextModelIntersection.idface == -1) return false; testParticle.nextModelIntersection.collisionPosition = testParticle.position + testParticle.direction*t; obst_dist += testParticle.direction.length()*t; From 9dfba91d6c3cafa5209f55eba2fadc1447fa6957 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 28 Oct 2021 16:57:05 +0200 Subject: [PATCH 33/56] Don't cast shadow rays with zero energy --- src/spps-agh/calculation_cores/NextEventEstimationCore.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp index 94acebb9..7d6cb6a0 100644 --- a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp +++ b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp @@ -336,6 +336,8 @@ void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t energy = BRDFs::SolveBRDFReflection(*materialInfo, faceNormal, receiver->position, shadowRay, particle.direction, configurationTool); } + if (energy == 0) break; + shadowRay.energie *= energy; //fast forward particle to receiver surrounding From 43cfdd9b7c7a4b7fe3b949a1829ceff58f57d123 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 28 Oct 2021 18:01:08 +0200 Subject: [PATCH 34/56] Fix rare case when the result is zeros only --- src/spps-agh/calculation_cores/MLTCore.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index c67dba81..95f6a8b0 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -76,6 +76,7 @@ void MLTCore::CalculateB() { b += seed.probability.totalProbability(wDC) / totalSeedNumber; } if (isnan(b)) std::cout << "!!!!!! b is NaN !!!!!!!!!" << std::endl; + if (b==0) std::cout << "!!!!!! b is zero, the result will be zeros only, use more seeds !!!!!!!!!" << std::endl; } void MLTCore::ChooseSeeds() { @@ -149,6 +150,9 @@ void MLTCore::RunMutation() void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) { CONF_PARTICULE_AGH shadowRay; + if (input_particle.weight == 0) { + input_particle.shadowRays.clear(); + } while (!input_particle.shadowRays.empty()) { @@ -156,7 +160,7 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) input_particle.shadowRays.pop_front(); //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; - if (selectSeeds) shadowRay.energie *= InitialSeedsSize / SeedsSize; // Compensate energy for seed selection + if (selectSeeds) shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection Run(shadowRay); //applicationTools.outputTool->SaveParticule(); } From 915da614b2b8c6c813f68d88d4258b672956d9e1 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 29 Oct 2021 00:36:54 +0200 Subject: [PATCH 35/56] Fix occasional crash on visability check - pt.2 --- src/spps-agh/CalculationCoreAGH.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spps-agh/CalculationCoreAGH.cpp b/src/spps-agh/CalculationCoreAGH.cpp index f218be76..46bd1687 100644 --- a/src/spps-agh/CalculationCoreAGH.cpp +++ b/src/spps-agh/CalculationCoreAGH.cpp @@ -201,6 +201,7 @@ bool CalculationCoreSPPS::VisabilityTest(const CONF_PARTICULE_AGH &configuration CONF_PARTICULE_AGH testParticle = configurationP; testParticle.nextModelIntersection.idface = GetTetraFaceCollision(testParticle, testParticle.direction, t); + if (testParticle.nextModelIntersection.idface == -1) return false; testParticle.nextModelIntersection.collisionPosition = testParticle.position + testParticle.direction*t; obst_dist = configurationP.direction.length()*t; From ddfd08ca532b3cdb5db63788538a903b873629aa Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 11 Nov 2021 18:03:46 +0100 Subject: [PATCH 36/56] Fix info message --- src/spps-agh/data_manager/core_configurationAGH.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index bfd6203a..19b9679f 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -222,7 +222,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Large step probability: " << largeStepProb << std::endl; //std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << // "Currently spec. ref. multiplier is disabled" << std::endl; - std::cout << "Refection penalty multiplier :" << reflPenaltyMultip << std::endl; + std::cout << "Reflection penalty multiplier: " << reflPenaltyMultip << std::endl; std::cout << "Seed selection: " << selectSeeds << std::endl; if (selectSeeds) std::cout << " Seed to sample rato:" << seedSampleRatio << std::endl; if (selectSeeds) std::cout << " Seed selection exponent:" << seedSelectionExponent << std::endl; From 2bb67e81437e7ee1cf6e85cc5d08be5c3ccd542d Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 11 Nov 2021 18:03:55 +0100 Subject: [PATCH 37/56] Increment version --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 04490b94..4a9a86c2 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -32,7 +32,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.3 (22.10.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.5 (11.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 0cdfc3c60eac1e8fc63f06ee504a45c96d4d9377 Mon Sep 17 00:00:00 2001 From: wbinek Date: Mon, 15 Nov 2021 13:24:16 +0100 Subject: [PATCH 38/56] Make GenerateShadowRays template to accept both vector and list argument --- .../NextEventEstimationCore.cpp | 10 +++++-- .../NextEventEstimationCore.h | 4 ++- src/spps-agh/sppsNeeAGHTypes.h | 2 +- src/spps-agh/tools/brdfreflection.cpp | 28 +++++++++---------- src/spps-agh/tools/brdfreflection.h | 10 +++---- 5 files changed, 30 insertions(+), 24 deletions(-) diff --git a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp index 7d6cb6a0..d6a5cf7a 100644 --- a/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp +++ b/src/spps-agh/calculation_cores/NextEventEstimationCore.cpp @@ -307,7 +307,8 @@ void NextEventEstimationCore::FreeParticleTranslation(CONF_PARTICULE_AGH &config configurationP.position += translationVector; } -void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo,const vec3& faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace* faceInfo, std::vector* probability) +template +void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo,const vec3& faceNormal,const double& deltaT,const double& distanceToTravel, C& shadowRays, const t_cFace* faceInfo, std::vector* probability) { //Calculate and cast shadow rays for each (t_Recepteur_P* receiver in configurationTool->recepteur_p_List) @@ -333,7 +334,7 @@ void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t energy = faceInfo->faceMaterial->customBrdf->getEnergy(curentFreq, faceNormal, particle.direction, newDirection)* solidAngle*0.66; } else { - energy = BRDFs::SolveBRDFReflection(*materialInfo, faceNormal, receiver->position, shadowRay, particle.direction, configurationTool); + energy = BRDFs::SolveBRDFReflection(*materialInfo, faceNormal, receiver->position, shadowRay.position, particle.direction, configurationTool); } if (energy == 0) break; @@ -362,4 +363,7 @@ void NextEventEstimationCore::GenerateShadowRays(CONF_PARTICULE_AGH& particle, t shadowRays.push_back(shadowRay); } } -} \ No newline at end of file +} + +template void NextEventEstimationCore::GenerateShadowRays>(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo, const vec3& faceNormal, const double& deltaT, const double& distanceToTravel, std::vector& shadowRays, const t_cFace* faceInfo, std::vector* probability); +template void NextEventEstimationCore::GenerateShadowRays>(CONF_PARTICULE_AGH& particle, t_Material_BFreq* materialInfo, const vec3& faceNormal, const double& deltaT, const double& distanceToTravel, std::list& shadowRays, const t_cFace* faceInfo, std::vector* probability); \ No newline at end of file diff --git a/src/spps-agh/calculation_cores/NextEventEstimationCore.h b/src/spps-agh/calculation_cores/NextEventEstimationCore.h index 576295e2..1ec5303c 100644 --- a/src/spps-agh/calculation_cores/NextEventEstimationCore.h +++ b/src/spps-agh/calculation_cores/NextEventEstimationCore.h @@ -13,7 +13,9 @@ class NextEventEstimationCore : public CalculationCoreSPPS protected: virtual void Movement(CONF_PARTICULE_AGH &configurationP) override; virtual void FreeParticleTranslation(CONF_PARTICULE_AGH &configurationP, const vec3 &translationVector) override; - void GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq * materialInfo, const vec3 & faceNormal,const double& deltaT,const double& distanceToTravel, std::list& shadowRays, const t_cFace * faceInfo, std::vector* probability = nullptr); + + template + void GenerateShadowRays(CONF_PARTICULE_AGH& particle, t_Material_BFreq * materialInfo, const vec3 & faceNormal,const double& deltaT,const double& distanceToTravel, C& shadowRays, const t_cFace * faceInfo, std::vector* probability = nullptr); }; #endif \ No newline at end of file diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 4a9a86c2..9de42e75 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -63,7 +63,7 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH std::vector refDir2; double sourceDir1; double sourceDir2; - std::list shadowRays; + std::vector shadowRays; double totalDistance = 0; double weight = 0; diff --git a/src/spps-agh/tools/brdfreflection.cpp b/src/spps-agh/tools/brdfreflection.cpp index ea11486e..6bbf3390 100644 --- a/src/spps-agh/tools/brdfreflection.cpp +++ b/src/spps-agh/tools/brdfreflection.cpp @@ -71,21 +71,21 @@ double BRDFs::calculatePNF(double d, int n) { } -double BRDFs::SolveBRDFReflection(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH &shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH *configurationTool) +double BRDFs::SolveBRDFReflection(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH *configurationTool) { switch (material.reflectionLaw) { case REFLECTION_LAW_SPECULAR: - return SolveSpecularBRDF(material, faceNormal, targetPoint, shadowRay, incomingDirection, configurationTool); + return SolveSpecularBRDF(material, faceNormal, targetPoint, position, incomingDirection, configurationTool); case REFLECTION_LAW_LAMBERT: - return SolveSpecularLambertBRDF(material, faceNormal, targetPoint, shadowRay, incomingDirection, configurationTool); + return SolveSpecularLambertBRDF(material, faceNormal, targetPoint, position, incomingDirection, configurationTool); case REFLECTION_LAW_PHONG: - return SolvePhongBRDF(material, faceNormal, targetPoint, shadowRay, incomingDirection, configurationTool); + return SolvePhongBRDF(material, faceNormal, targetPoint, position, incomingDirection, configurationTool); default: - return SolveSpecularLambertBRDF(material, targetPoint, faceNormal, shadowRay, incomingDirection, configurationTool); + return SolveSpecularLambertBRDF(material, targetPoint, faceNormal, position, incomingDirection, configurationTool); } } @@ -95,14 +95,14 @@ vec3 BRDFs::SolveSpecularReflection(const vec3 &incomingDirection, const vec3 &f return retVal / retVal.length(); } -double BRDFs::SolveSpecularBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH& shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) +double BRDFs::SolveSpecularBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) { vec3 specular = SolveSpecularReflection(incomingDirection, faceNormal); double specularEnergy = 0; l_decimal mu1, mu2; float receiverRadius = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_RAYON_RECEPTEURP); - if (RaySphereIntersection(shadowRay.position, shadowRay.position + specular * 1000, targetPoint, receiverRadius, &mu1, &mu2)) + if (RaySphereIntersection(position, position + specular * 1000, targetPoint, receiverRadius, &mu1, &mu2)) { specularEnergy = (1 - material.diffusion) * (abs(mu2 - mu1) * specular.length() * 1000) / (2 * receiverRadius); } @@ -110,15 +110,15 @@ double BRDFs::SolveSpecularBRDF(const t_Material_BFreq& material, const vec3& fa return specularEnergy; } -double BRDFs::SolveSpecularLambertBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH& shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) +double BRDFs::SolveSpecularLambertBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) { vec3 specular = SolveSpecularReflection(incomingDirection, faceNormal); - vec3 toReceiver = targetPoint - shadowRay.position; + vec3 toReceiver = targetPoint - position; double specularEnergy = 0, lambertEnergy; l_decimal mu1, mu2; float receiverRadius = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_RAYON_RECEPTEURP); - if (RaySphereIntersection(shadowRay.position, shadowRay.position + specular * 1000, targetPoint, receiverRadius, &mu1, &mu2)) + if (RaySphereIntersection(position, position + specular * 1000, targetPoint, receiverRadius, &mu1, &mu2)) { specularEnergy = (1 - material.diffusion) * (abs(mu2 - mu1) * specular.length()*1000)/(2*receiverRadius); } @@ -135,10 +135,10 @@ double BRDFs::SolveSpecularLambertBRDF(const t_Material_BFreq& material, const v return specularEnergy + lambertEnergy; } -double BRDFs::SolvePhongBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH& shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) +double BRDFs::SolvePhongBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool) { vec3 specular = SolveSpecularReflection(incomingDirection, faceNormal); - vec3 toReceiver = targetPoint - shadowRay.position; + vec3 toReceiver = targetPoint - position; Matrix3 rotMatrix; double energyFactor = 0, solidAngle; float receiverRadius = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_RAYON_RECEPTEURP); @@ -150,12 +150,12 @@ double BRDFs::SolvePhongBRDF(const t_Material_BFreq& material, const vec3& faceN // consider taking n into accont - small n means slowly changing Phong function. solidAngle = (M_PI*receiverRadius*receiverRadius) / (toReceiver.length()*toReceiver.length()); - evaluatePhongAtPoint(n, material.diffusion, solidAngle, targetPoint, shadowRay.position, faceNormal, specular, energyFactor); + evaluatePhongAtPoint(n, material.diffusion, solidAngle, targetPoint, position, faceNormal, specular, energyFactor); return energyFactor*0.66; } -void BRDFs::evaluatePhongAtPoint(int n, float diffusion, double solidAngle, vec3 target, vec3 position, vec3 faceNormal, vec3 specular, double &result) +void BRDFs::evaluatePhongAtPoint(const int n, const float diffusion, const double solidAngle, const vec3& target, const vec3& position, const vec3& faceNormal, const vec3& specular, double &result) { vec3 toSubSurf = target - position; diff --git a/src/spps-agh/tools/brdfreflection.h b/src/spps-agh/tools/brdfreflection.h index 563415d9..8bf4749b 100644 --- a/src/spps-agh/tools/brdfreflection.h +++ b/src/spps-agh/tools/brdfreflection.h @@ -28,7 +28,7 @@ class BRDFs * @param incomingDirection Direction of incoming ray * @param configurationTool Simulation configuration */ - static double SolveBRDFReflection(const t_Material_BFreq& material, const vec3& faceNormal,const vec3& targetPoint, const CONF_PARTICULE_AGH &shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH *configurationTool); + static double SolveBRDFReflection(const t_Material_BFreq& material, const vec3& faceNormal,const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH *configurationTool); /** * Calculates specular reflection direction @@ -38,7 +38,7 @@ class BRDFs */ static vec3 SolveSpecularReflection(const vec3& incomingDirection, const vec3& faceNormal); - static double SolveSpecularBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH & shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH * configurationTool); + static double SolveSpecularBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH * configurationTool); /** * Solves reflecion energy using Specular+Lambert BRDF @@ -49,7 +49,7 @@ class BRDFs * @param incomingDirection Direction of incoming ray * @param configurationTool Simulation configuration */ - static double SolveSpecularLambertBRDF(const t_Material_BFreq& material,const vec3& faceNormal,const vec3& targetPoint, const CONF_PARTICULE_AGH& shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool); + static double SolveSpecularLambertBRDF(const t_Material_BFreq& material,const vec3& faceNormal,const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool); /** * Solves reflecion energy using Phong BRDF @@ -61,8 +61,8 @@ class BRDFs * @param incomingDirection Direction of incoming ray * @param configurationTool Simulation configuration */ - static double SolvePhongBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const CONF_PARTICULE_AGH& shadowRay, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool); - static void evaluatePhongAtPoint(int n, float diffusion, double solidAngle, vec3 target, vec3 position, vec3 faceNormal, vec3 specular, double& result); + static double SolvePhongBRDF(const t_Material_BFreq& material, const vec3& faceNormal, const vec3& targetPoint, const vec3& position, const vec3& incomingDirection, Core_ConfigurationAGH* configurationTool); + static void evaluatePhongAtPoint(const int n, const float diffusion, const double solidAngle, const vec3& target, const vec3& position, const vec3& faceNormal, const vec3& specular, double& result); static double calculatePhongNormalizationFactor(const vec3 & faceNormal, const vec3 & specularDirection, const int & n); }; \ No newline at end of file From 3689a12962791982f2ba93266425e6ca72d0d1a1 Mon Sep 17 00:00:00 2001 From: wbinek Date: Mon, 15 Nov 2021 13:24:59 +0100 Subject: [PATCH 39/56] Add low memory mode for MLT algorithm --- .../tree_core/e_core_spps_aghcore_advanced.h | 1 + src/spps-agh/calculation_cores/MLTCore.cpp | 12 ++++++++++-- src/spps-agh/calculation_cores/MLTCore.h | 2 ++ src/spps-agh/data_manager/core_configurationAGH.cpp | 6 ++++++ src/spps-agh/data_manager/core_configurationAGH.h | 1 + 5 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index b2b68909..c6d2bf3f 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -293,6 +293,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); confCore->AppendPropertyDecimal("seedSelectionExponent", wxTRANSLATE("Tune - Seed selection exponent"), 0.5, false, 2, true, true, 1, 0.01, true); + confCore->AppendPropertyBool("lowMemoryMode", wxTRANSLATE("Tune - Low memory mode"), true, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 95f6a8b0..d3eb8eda 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -30,6 +30,8 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT genS1 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S1); genS2 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S2); + + lowMemoryMode = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_LOW_MEMORY_MODE); }; bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) @@ -46,6 +48,8 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) inputParticle.weight = 0; RunInitialSeed(inputParticle); + + if (lowMemoryMode) inputParticle.shadowRays.clear(); InitialSeeds.push_back(inputParticle); } else @@ -143,6 +147,10 @@ void MLTCore::RunMutation() { CONF_PARTICULE_MLT inputParticle = Seeds.front(); Seeds.pop_front(); + if (lowMemoryMode) { + inputParticle.weight = 0; + Follow(inputParticle); + } inputParticle.isInitialSeed = false; DoMutations(inputParticle, (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size())); } @@ -156,8 +164,8 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) while (!input_particle.shadowRays.empty()) { - shadowRay = input_particle.shadowRays.front(); - input_particle.shadowRays.pop_front(); + shadowRay = input_particle.shadowRays.back(); + input_particle.shadowRays.pop_back(); //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; if (selectSeeds) shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index b10b6229..55c9257d 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -35,6 +35,8 @@ class MLTCore : public NextEventEstimationCore bool mutateSource = true; bool weightDecayCompensation = true; + bool lowMemoryMode = true; + int genS1 = 1024; int genS2 = 64; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 19b9679f..c7bc7e01 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -218,6 +218,8 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int mutationGeneratorS1 = advancedNode->GetProperty("mutationGeneratorS1").ToInt(); int mutationGeneratorS2 = advancedNode->GetProperty("mutationGeneratorS2").ToInt(); + int lowMemoryMode = advancedNode->GetProperty("lowMemoryMode").ToInt(); + std::cout << "Mutation number: " << mutationNumber << std::endl; std::cout << "Large step probability: " << largeStepProb << std::endl; //std::cout << "Specular reflection probability multiplier: " << specularReflProbabilityMulti << std::endl << @@ -238,6 +240,8 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Mutation generator S1: " << mutationGeneratorS1 << std::endl; std::cout << "Mutation generator S2: " << mutationGeneratorS2 << std::endl; + + std::cout << "Low memory mode: " << lowMemoryMode << std::endl; if (useWeightDecayCompensation && reflPenaltyMultip > 0.95) std::cout << "WARNING! Using Weight Decay Compensation with high relfection penalty multiplier may cause artefacts!" << std::endl; @@ -261,6 +265,8 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S1, mutationGeneratorS1); SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S2, mutationGeneratorS2); + + SetConfigInformation(IPROP_MLT_LOW_MEMORY_MODE, lowMemoryMode); } } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index 39c05c78..f5ec9fd0 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -56,6 +56,7 @@ class Core_ConfigurationAGH : public Core_Configuration IPROP_MLT_COST_AIR_ABSORPTION, IPROP_MLT_MUTATION_GENERATOR_S1, IPROP_MLT_MUTATION_GENERATOR_S2, + IPROP_MLT_LOW_MEMORY_MODE, }; /** * Initialisation des paramètres du coeur de calcul à partir d'un fichier XML From 4c93b5e6e4adf5b580780ed4e776cdbb512482ee Mon Sep 17 00:00:00 2001 From: wbinek Date: Mon, 15 Nov 2021 13:27:50 +0100 Subject: [PATCH 40/56] Increment SPPS-AGH version --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 9de42e75..0147a4fe 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -32,7 +32,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.5 (11.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.6 (15.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 40b6eb24dcccc1c2b2ecea166fda68c40689f8b9 Mon Sep 17 00:00:00 2001 From: wbinek Date: Mon, 15 Nov 2021 15:19:25 +0100 Subject: [PATCH 41/56] Add more memory optimizations, increment SPPS-AGH version --- src/spps-agh/calculation_cores/MLTCore.cpp | 15 +++++++++------ src/spps-agh/sppsNeeAGHTypes.h | 4 ++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index d3eb8eda..11d46b27 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -65,7 +65,7 @@ void MLTCore::CalculatewDC() { { int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); wDC = 0; - for (CONF_PARTICULE_MLT seed : InitialSeeds) { + for (const CONF_PARTICULE_MLT& seed : InitialSeeds) { double currentProbabilityRaw = seed.probability.currentProbability / pow(reflectionPenaltyMultip, seed.reflectionOrder); double wDCIncrement = pow(currentProbabilityRaw, (1. / (seed.reflectionsNum + 1))) / totalSeedNumber; wDC += wDCIncrement; @@ -76,7 +76,7 @@ void MLTCore::CalculatewDC() { void MLTCore::CalculateB() { int totalSeedNumber = (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size()); - for (CONF_PARTICULE_MLT seed : InitialSeeds) { + for (const CONF_PARTICULE_MLT& seed : InitialSeeds) { b += seed.probability.totalProbability(wDC) / totalSeedNumber; } if (isnan(b)) std::cout << "!!!!!! b is NaN !!!!!!!!!" << std::endl; @@ -88,7 +88,7 @@ void MLTCore::ChooseSeeds() { // If seed chosing disabled if (!selectSeeds) { - Seeds = InitialSeeds; + Seeds = std::move(InitialSeeds); std::cout << "Seed selection disabled, copied all initial seeds" << std::endl; return; } @@ -101,11 +101,12 @@ void MLTCore::ChooseSeeds() { int individualRays = 0; bool accounted = false; - for (CONF_PARTICULE_MLT seed : InitialSeeds) { + for (const CONF_PARTICULE_MLT& seed : InitialSeeds) { total_score += pow(seed.probability.totalProbability(wDC), seedSelectionExponent); } - for (CONF_PARTICULE_MLT seed : InitialSeeds) { + while (!InitialSeeds.empty()) { + CONF_PARTICULE_MLT seed = InitialSeeds.back(); accounted = false; float currentScore = pow(seed.probability.totalProbability(wDC), seedSelectionExponent) / total_score * (N_target); if (currentScore > 1) { @@ -123,6 +124,8 @@ void MLTCore::ChooseSeeds() { Seeds.push_back(seed); if (!accounted) individualRays++; } + // Release memory of the seed on the go + InitialSeeds.pop_back(); } SeedsSize = Seeds.size(); float avgRayRepeat = (float)SeedsSize / individualRays; @@ -135,7 +138,7 @@ void MLTCore::ChooseSeeds() { "Consider lower sample to seed ratio" << std::endl; } // Release memory taken by initial seeds, as they are no more reqired - InitialSeeds.clear(); + // InitialSeeds.clear(); } bool MLTCore::SeedIsEmpty() diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 0147a4fe..509c540f 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -32,7 +32,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.6 (15.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.7 (15.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION @@ -80,7 +80,7 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH std::vector partialProbability; //probability of generating shadow ray std::vector srRelfectionOrder; //Used for weight decay compensation - double totalProbability(double wDC) //sum of shadow rays probability - total ray contribution + double totalProbability(double wDC) const //sum of shadow rays probability - total ray contribution { double totalProb = 0; for (int i = 0; i < partialProbability.size();i++) { From 949e263ad6340264f1c2f80e27a3ee8bf3de9fdc Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 18 Nov 2021 18:49:05 +0100 Subject: [PATCH 42/56] MLT cast all initial rays (quality improvement), increment SPPS-AGH version --- src/spps-agh/calculation_cores/MLTCore.cpp | 54 +++++++++++----------- src/spps-agh/calculation_cores/MLTCore.h | 6 +-- src/spps-agh/sppsNeeAGHTypes.h | 12 +---- 3 files changed, 32 insertions(+), 40 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 11d46b27..d64422a1 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -31,24 +31,26 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT genS1 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S1); genS2 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S2); - lowMemoryMode = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_LOW_MEMORY_MODE); + lowMemoryMode = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_LOW_MEMORY_MODE); + + InitialSeedsSize = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION); + SeedsSize = selectSeeds ? (int)(InitialSeedsSize * seedSampleRatio) : InitialSeedsSize; }; bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) - -//TODO: Mutacja kierunku w którym generowany jest promieñ! - -//TODO: Przenieœæ initial seed przed mutacje, zapisaæ b, zastanowiæ siê nad stanem pocz¹tkowym! - DONE -//TODO: mno¿enie razy prawdopodobieñstwo SR w danym kierunku - powinno znacznie pomóc - DONE { if (!configurationP.isShadowRay) { CONF_PARTICULE_MLT inputParticle(configurationP); inputParticle.SetInitialState(inputParticle); - inputParticle.weight = 0; - RunInitialSeed(inputParticle); + // The ray will be casted as it comes from random distribution, the weight is assigned so that it caries right energy + // ((double)InitialSeedsSize/SeedsSize) is essentialy to cancel with weight assigment in CastShadowRay function + inputParticle.weight = 1./ (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize/SeedsSize)); + RunInitialSeed(inputParticle); + if (rejectionSampling) + CastShadowRays(inputParticle, lowMemoryMode); if (lowMemoryMode) inputParticle.shadowRays.clear(); InitialSeeds.push_back(inputParticle); } @@ -150,31 +152,29 @@ void MLTCore::RunMutation() { CONF_PARTICULE_MLT inputParticle = Seeds.front(); Seeds.pop_front(); - if (lowMemoryMode) { - inputParticle.weight = 0; - Follow(inputParticle); - } + inputParticle.weight = 0; + if (mutation_number == 0) return; + if (lowMemoryMode) Follow(inputParticle); + inputParticle.isInitialSeed = false; DoMutations(inputParticle, (*configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION) * configurationTool->srcList.size())); } -void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle) +void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle, bool delete_sr = true) { - CONF_PARTICULE_AGH shadowRay; if (input_particle.weight == 0) { input_particle.shadowRays.clear(); } - while (!input_particle.shadowRays.empty()) + for(CONF_PARTICULE_AGH& shadowRay : input_particle.shadowRays) { - shadowRay = input_particle.shadowRays.back(); - input_particle.shadowRays.pop_back(); //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; - if (selectSeeds) shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection + shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection, if no seed selection the sizes are equal Run(shadowRay); //applicationTools.outputTool->SaveParticule(); } + if (delete_sr) input_particle.shadowRays.clear(); } void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNumber) @@ -183,9 +183,11 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu double ratio; int mutation_counter = 0; - if (mutation_number == 0) - if (selectSeeds) inputParticle.weight = b / inputParticle.probability.totalProbability(wDC); //poprawne dla próbek po wyborze z puli - else inputParticle.weight = 1; //poprawne dla losowych probek + // We assure that thatu mutation number is greater than 0 in RunMutation function, further logic changed when seed casting for initial seeds was added, to be removed in future + //if (mutation_number == 0) + // return; + // if (selectSeeds) inputParticle.weight = b / inputParticle.probability.totalProbability(wDC); //poprawne dla próbek po wyborze z puli + // else inputParticle.weight = 1; //poprawne dla losowych probek while (mutation_counter < mutation_number) { @@ -195,20 +197,18 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu ratio = std::min(mutatedParticle.probability.totalProbability(wDC) / inputParticle.probability.totalProbability(wDC), 1.); if (isnan(ratio)) [[unlikely]] std::cout << "Ratio is NaN, may cause serious erros" << std::endl; - //inputParticle.weight += (1 - ratio) / ((inputParticle.totalProbability / b + pLarge)*(mutation_number)); - //mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / ((mutatedParticle.totalProbability / b + pLarge)*(mutation_number)); //// The if block calculate sample weights for diferent cases // Treat sample that is not large step if (rejectionSampling){ // Rejection sampling - inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1 - ratio)) / (mutation_number * (1 + pLarge)) ; - mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio) + mutatedParticle.largeStep) / (mutation_number * (1 + pLarge)); + inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio)) / (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize / SeedsSize)) ; + mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio) + mutatedParticle.largeStep) / (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize / SeedsSize)); } else{ // No rejection sampling - inputParticle.weight = b / inputParticle.probability.totalProbability(wDC) / (1 + pLarge); - mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability(wDC) + mutatedParticle.largeStep) / (1 + pLarge); + inputParticle.weight = b / inputParticle.probability.totalProbability(wDC) / (1. + pLarge); + mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability(wDC) + mutatedParticle.largeStep) / (1. + pLarge); } //// Make particle seletion, update and SR casting diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 55c9257d..fde3f848 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -44,12 +44,12 @@ class MLTCore : public NextEventEstimationCore std::list InitialSeeds; std::list Seeds; - int InitialSeedsSize; - int SeedsSize; + int InitialSeedsSize=0; + int SeedsSize=0; void Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_MLT& mutetedParticle) const; bool MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rndDir1, double rndDir2, double rndTravelProbability, double rndMatAbsorbtion, float deltaT, float distanceToTravel); - void CastShadowRays(CONF_PARTICULE_MLT& input_particle); + void CastShadowRays(CONF_PARTICULE_MLT& input_particle, bool delete_sr); void DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNumber); void Movement(CONF_PARTICULE_MLT &configurationP); void Follow(CONF_PARTICULE_MLT& propagationParticle); diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 509c540f..57fcb8a2 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -11,15 +11,7 @@ * ---------------------------------------------------------------------------------------------------------------*/ /* --------------------------------------------------------------------------------------------------------- -*changelog 10.09.2016: -* - Fixed: Next Event Estimation Code - passed evaluation -* - New: New calculation code - MLT (Kelemen style) - intial version -* -* - TODO: ALL - Full integration with source directivity options -* - TODO: MLT - Fix for models with more than one receiver - wrong weight calculation -* - TODO: MLT - Configuration via I-SIMPA interface, not hardcoded values -* - TODO: MLT - Initial ray direction mutations -* - TODO: MLT - Compliance with all diffusion models (currently only Lambert + Specular) + * ---------------------------------------------------------------------------------------------------------------*/ //#include "coreTypes.h" @@ -32,7 +24,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.4.7 (15.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.0beta (18.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From e058f7c98051a910bab540101385a36ff73b3f1b Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 18 Nov 2021 18:49:37 +0100 Subject: [PATCH 43/56] Disable most MLT settings if no mutations is 0 --- .../tree_core/e_core_spps_aghcore_advanced.h | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index c6d2bf3f..81a0504a 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -258,6 +258,22 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element } else if (filsInfo.libelleElement == "mutationNumber") { this->SetReadOnlyConfig("largeStepProb", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("selectSeeds", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("seedSampleRatio", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("seedSelectionExponent", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("useWeightDecayCompensation", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("mutateSource", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("reflectionPenalty", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("costMatAbsorption", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("costMatDiffusion", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("costShadowRays", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("costReflDir", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("costAirAb", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("useRejectionSampling", !(bool)this->GetIntegerConfig("mutationNumber")); + this->SetReadOnlyConfig("lowMemoryMode", !(bool)this->GetIntegerConfig("mutationNumber")); + + if (this->GetIntegerConfig("mutationNumber") == 0) this->UpdateBoolConfig("useRejectionSampling", true); + if (this->GetIntegerConfig("mutationNumber") == 0) this->UpdateBoolConfig("lowMemoryMode", true); } Element::Modified(eModif); @@ -278,9 +294,9 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element } void InitBaseMLTProperties(Element* confCore) { - confCore->AppendPropertyDecimal("largeStepProb", wxTRANSLATE("Tune - Large step probability"), 0.5, false, 2, true, true, 1, 0, true); + confCore->AppendPropertyDecimal("largeStepProb", wxTRANSLATE("Tune - Large step probability"), 0.2, false, 2, true, true, 1, 0, true); //confCore->AppendPropertyDecimal("specularReflProbabilityMulti", wxTRANSLATE("Tune - Specular reflection probability multiplier"), 0.5, false, 3, true, true, 10, 0.01, true); - confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Sampling - Mutation number"), 25, true, false, true, 0, 0); + confCore->AppendPropertyInteger("mutationNumber", wxTRANSLATE("Sampling - Mutation number"), 10, true, false, true, 0, 0); } void InitAdditionalMLTProperties(Element* confCore) { confCore->AppendPropertyBool("useRejectionSampling", wxTRANSLATE("Use - Rejection sampling"), true, true); @@ -296,9 +312,9 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyBool("lowMemoryMode", wxTRANSLATE("Tune - Low memory mode"), true, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); - confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); + confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, false); confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); - confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), true, true); + confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), true, false); confCore->AppendPropertyBool("costAirAbs", wxTRANSLATE("Cost - use air absorption"), false, true); } }; From ae5175f5867bd92402f5132fd2dc3fb430b2fdce Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 18 Nov 2021 22:13:41 +0100 Subject: [PATCH 44/56] MLT auto enable low memory mode and rejection sampling for 0 mutations --- src/spps-agh/calculation_cores/MLTCore.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index d64422a1..ca30e0d4 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -35,6 +35,12 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT InitialSeedsSize = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION); SeedsSize = selectSeeds ? (int)(InitialSeedsSize * seedSampleRatio) : InitialSeedsSize; + + if (mutation_number == 0) { + lowMemoryMode = true; + rejectionSampling = true; + std::cout << "Zero mutations - enabled low memory mode and rejection sampling, can't simulate otherwise" << std::endl; + } }; bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) From fe7775d56c492e0e69d6e5a940760f8b791e79ff Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 23 Nov 2021 23:30:21 +0100 Subject: [PATCH 45/56] Add single value mutation probability parameter --- .../tree_core/e_core_spps_aghcore_advanced.h | 1 + src/spps-agh/calculation_cores/MLTCore.cpp | 19 ++++++++++++------- src/spps-agh/calculation_cores/MLTCore.h | 3 ++- .../data_manager/core_configurationAGH.cpp | 3 +++ .../data_manager/core_configurationAGH.h | 1 + src/spps-agh/tools/dotreflectionAGH.h | 2 +- 6 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 81a0504a..2c048211 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -308,6 +308,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.8, false, 2, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); + confCore->AppendPropertyDecimal("mutationProbability", wxTRANSLATE("Tune - Mutation probability"), 1, false, 3, true, true, 1, 0.1, true); confCore->AppendPropertyDecimal("seedSelectionExponent", wxTRANSLATE("Tune - Seed selection exponent"), 0.5, false, 2, true, true, 1, 0.01, true); confCore->AppendPropertyBool("lowMemoryMode", wxTRANSLATE("Tune - Low memory mode"), true, true); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index ca30e0d4..26095436 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -30,6 +30,7 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT genS1 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S1); genS2 = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_MUTATION_GENERATOR_S2); + mutationProbability = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATION_PROBABILITY); lowMemoryMode = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_LOW_MEMORY_MODE); @@ -307,14 +308,14 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML mutetedParticle.sourceDir2 = GetRandValue(); }else{ - MutationGenerator(mutetedParticle.sourceDir1); - MutationGenerator(mutetedParticle.sourceDir2); + MutationGenerator(mutetedParticle.sourceDir1,0); + MutationGenerator(mutetedParticle.sourceDir2,0); for (int i = 0; i < input_particle.reflectionsNum; i++) { - MutationGenerator(mutetedParticle.refDir1[i]); - MutationGenerator(mutetedParticle.refDir2[i]); - MutationGenerator(mutetedParticle.travelProbability[i]); - MutationGenerator(mutetedParticle.matAbsorbtions[i]); + MutationGenerator(mutetedParticle.refDir1[i],i); + MutationGenerator(mutetedParticle.refDir2[i],i); + MutationGenerator(mutetedParticle.travelProbability[i], i); + MutationGenerator(mutetedParticle.matAbsorbtions[i], i); } } } @@ -550,8 +551,12 @@ bool MLTCore::GetNextCollision(CONF_PARTICULE_MLT& configurationP, double& dista } -void MLTCore::MutationGenerator(double& value) const +void MLTCore::MutationGenerator(double& value, const int& reflOrd) const { + float mut_prob = pow(mutationProbability, reflOrd+1); + if (GetRandValue() > mut_prob) { + return; + } float s1 = 1. / genS1, s2 = 1. / genS2; float dv = s2*exp(-log(s2 / s1)*GetRandValue()); if (GetRandValue() < 0.5) { diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index fde3f848..bff2512b 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -39,6 +39,7 @@ class MLTCore : public NextEventEstimationCore int genS1 = 1024; int genS2 = 64; + float mutationProbability = 1.f; protected: std::list InitialSeeds; @@ -55,7 +56,7 @@ class MLTCore : public NextEventEstimationCore void Follow(CONF_PARTICULE_MLT& propagationParticle); void FreeParticleTranslation(CONF_PARTICULE_AGH &configurationP, const vec3 &translationVector) override; bool GetNextCollision(CONF_PARTICULE_MLT &configurationP, double &distance); - void MutationGenerator(double& value) const; + void MutationGenerator(double& value, const int& reflOrd) const; }; #endif \ No newline at end of file diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index c7bc7e01..5ed1ce2b 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -217,6 +217,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int mutationGeneratorS1 = advancedNode->GetProperty("mutationGeneratorS1").ToInt(); int mutationGeneratorS2 = advancedNode->GetProperty("mutationGeneratorS2").ToInt(); + float mutationProbability = advancedNode->GetProperty("mutationProbability").ToFloat(); int lowMemoryMode = advancedNode->GetProperty("lowMemoryMode").ToInt(); @@ -240,6 +241,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Mutation generator S1: " << mutationGeneratorS1 << std::endl; std::cout << "Mutation generator S2: " << mutationGeneratorS2 << std::endl; + std::cout << "Mutation probability: " << mutationProbability << std::endl; std::cout << "Low memory mode: " << lowMemoryMode << std::endl; @@ -265,6 +267,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S1, mutationGeneratorS1); SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S2, mutationGeneratorS2); + SetConfigInformation(FPROP_MLT_MUTATION_PROBABILITY, mutationProbability); SetConfigInformation(IPROP_MLT_LOW_MEMORY_MODE, lowMemoryMode); } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index f5ec9fd0..b0048c94 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -33,6 +33,7 @@ class Core_ConfigurationAGH : public Core_Configuration FPROP_MLT_SEED_SAMPLE_RATIO, FPROP_MLT_REFL_PENALTY_MULTIP, FPROP_MLT_SEED_SELECTION_EXPONENT, + FPROP_MLT_MUTATION_PROBABILITY, }; enum NEW_IPROP diff --git a/src/spps-agh/tools/dotreflectionAGH.h b/src/spps-agh/tools/dotreflectionAGH.h index 0b2381da..81d6b309 100644 --- a/src/spps-agh/tools/dotreflectionAGH.h +++ b/src/spps-agh/tools/dotreflectionAGH.h @@ -136,7 +136,7 @@ class ReflectionLawsMLT : public ReflectionLawsAGH static vec3 PhongSpecularPart(vec3 &vectorDirection, vec3 &faceNormal, CONF_PARTICULE_AGH& particuleInfo, t_Material_BFreq &material, double &rnd1, double &rnd2, double &probability) { - float n = powf(10, -1.7234170470604733 * material.diffusion + 2.6245274102195886); + int n = (int) powf(10, -1.7234170470604733 * material.diffusion + 2.6245274102195886); vec3 specularDir = SpecularReflection(vectorDirection, faceNormal); Matrix3 rotationMatrix; From 6214444b9db1aadb5ac5651f2634ec20e1eb95d4 Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 23 Nov 2021 23:31:01 +0100 Subject: [PATCH 46/56] Increment spps-agh version to 0.5.1beta --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 57fcb8a2..547cd4b6 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -24,7 +24,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.0beta (18.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.1beta (23.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 9dc1b3cf28b9ae2c0633a4968510add8e0b859eb Mon Sep 17 00:00:00 2001 From: wbinek Date: Tue, 23 Nov 2021 23:56:42 +0100 Subject: [PATCH 47/56] Fix export to core settings for MLT --- .../data_manager/tree_core/e_core_spps_aghcore_advanced.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 2c048211..39c4cfe0 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -313,9 +313,9 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyBool("lowMemoryMode", wxTRANSLATE("Tune - Low memory mode"), true, true); confCore->AppendPropertyBool("costMatAbsorption", wxTRANSLATE("Cost - use material absorption"), true, true); - confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, false); + confCore->AppendPropertyBool("costMatDiffusion", wxTRANSLATE("Cost - use material diffusion"), false, true); confCore->AppendPropertyBool("costShadowRays", wxTRANSLATE("Cost - use shadow rays"), true, true); - confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), true, false); + confCore->AppendPropertyBool("costReflDir", wxTRANSLATE("Cost - use reflection direction"), false, true); confCore->AppendPropertyBool("costAirAbs", wxTRANSLATE("Cost - use air absorption"), false, true); } }; From c52b742b23520b14cb813c9fd08019f09d891d63 Mon Sep 17 00:00:00 2001 From: wbinek Date: Wed, 24 Nov 2021 00:22:44 +0100 Subject: [PATCH 48/56] Change mutationProbability, reflection order dependance --- src/spps-agh/calculation_cores/MLTCore.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 26095436..a46df948 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -312,10 +312,10 @@ void MLTCore::Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_ML MutationGenerator(mutetedParticle.sourceDir2,0); for (int i = 0; i < input_particle.reflectionsNum; i++) { - MutationGenerator(mutetedParticle.refDir1[i],i); - MutationGenerator(mutetedParticle.refDir2[i],i); - MutationGenerator(mutetedParticle.travelProbability[i], i); - MutationGenerator(mutetedParticle.matAbsorbtions[i], i); + MutationGenerator(mutetedParticle.refDir1[i],i+1); + MutationGenerator(mutetedParticle.refDir2[i],i+1); + MutationGenerator(mutetedParticle.travelProbability[i], i+1); + MutationGenerator(mutetedParticle.matAbsorbtions[i], i+1); } } } @@ -553,7 +553,7 @@ bool MLTCore::GetNextCollision(CONF_PARTICULE_MLT& configurationP, double& dista void MLTCore::MutationGenerator(double& value, const int& reflOrd) const { - float mut_prob = pow(mutationProbability, reflOrd+1); + float mut_prob = pow(mutationProbability, reflOrd); if (GetRandValue() > mut_prob) { return; } From 845891047636d1b0f306572f46021e4efc972775 Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 7 Jan 2022 17:30:26 +0100 Subject: [PATCH 49/56] Fix occasional NaNs in LF and LFC calculation --- src/lib_interface/Core/mathlib.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/lib_interface/Core/mathlib.h b/src/lib_interface/Core/mathlib.h index 62bee127..f5d0944e 100644 --- a/src/lib_interface/Core/mathlib.h +++ b/src/lib_interface/Core/mathlib.h @@ -179,7 +179,10 @@ class base_vec3 { // this < _v bool compare(const base_vec3 &_v,base_t epsi=EPSILON) { return (fabs(this->x - _v.x) < epsi && fabs(this->y - _v.y) < epsi && fabs(this->z - _v.z) < epsi); } base_t angle(const base_vec3 &v) const { // retourne l'angle en radians entre *this et v - base_t angle = acos(this->dot(v)/(this->length()*v.length())); + base_t dot_prod = this->dot(v) / (this->length() * v.length()); + if (dot_prod > 1) dot_prod = 1; + if (dot_prod < -1) dot_prod = -1; + base_t angle = acos(dot_prod); if(angle < EPSILON) return 0; return angle; } From b89dd570e3db0a1453e975a38d3a9b807f61b1bd Mon Sep 17 00:00:00 2001 From: wbinek Date: Fri, 7 Jan 2022 17:31:34 +0100 Subject: [PATCH 50/56] Add options to add 10 and 50 calculations to job_tool --- .../SystemScript/job_tool/__init__.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/currentRelease/SystemScript/job_tool/__init__.py b/currentRelease/SystemScript/job_tool/__init__.py index 2e52f090..ae3aee7a 100644 --- a/currentRelease/SystemScript/job_tool/__init__.py +++ b/currentRelease/SystemScript/job_tool/__init__.py @@ -3,7 +3,7 @@ import operator import uilocale import os -import collections +from collections.abc import Sequence _=uilocale.InstallUiModule(ui.application.getapplicationpath()["systemscript"]+"job_tool"+os.sep,ui.application.getlocale()) @@ -99,7 +99,7 @@ class JobManager(object): """ def appendjob(self,jobname,identifier): #if operator.isSequenceType(identifier): - if isinstance(identifier,collections.Sequence): + if isinstance(identifier,Sequence): self.joblst.append(job_types[jobname](identifier)) else: self.joblst.append(job_types[jobname](jobname,identifier)) @@ -118,6 +118,16 @@ def AddRun(self,idel): self.AddOpenProjectIfChanges() self.appendjob("RUN",idel) self.appendjob("SAVE_PROJECT",idel) + def AddRun10(self,idel): + self.AddOpenProjectIfChanges() + for i in range(10): + self.appendjob("RUN",idel) + self.appendjob("SAVE_PROJECT",idel) + def AddRun50(self,idel): + self.AddOpenProjectIfChanges() + for i in range(50): + self.appendjob("RUN",idel) + self.appendjob("SAVE_PROJECT",idel) def ClearJobLst(self,idel): self.joblst=[] self.lastprojectpath="" @@ -134,12 +144,16 @@ class manager: def __init__(self): self.jmanager=JobManager() self.addcalculation_id=ui.application.register_event(self.jmanager.AddRun) + self.addcalculation10_id=ui.application.register_event(self.jmanager.AddRun10) + self.addcalculation50_id=ui.application.register_event(self.jmanager.AddRun50) self.clear_jobs_id=ui.application.register_event(self.jmanager.ClearJobLst) self.execute_jobs_id=ui.application.register_event(self.jmanager.ExecJobs) self.print_jobs_id=ui.application.register_event(self.jmanager.PrintJobLst) def getmenu(self,typeel,idel,menu): submenu=[] submenu.append((_(u"Add calculation in the job list"),self.addcalculation_id)) + submenu.append((_(u"Add 10 calculations to the job list"),self.addcalculation10_id)) + submenu.append((_(u"Add 50 calculations to the job list"),self.addcalculation50_id)) submenu.append(()) submenu.append((_(u"Clear job list"),[(_(u"Confirmation"),self.clear_jobs_id)])) submenu.append((_(u"Run job list"),self.execute_jobs_id)) From bedb75885b10d84183c7b299027d17629cf5f3e6 Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 8 Jan 2022 11:02:31 +0100 Subject: [PATCH 51/56] Add mutated particle weight parameter to MLT algorithm --- .../data_manager/tree_core/e_core_spps_aghcore_advanced.h | 1 + src/spps-agh/calculation_cores/MLTCore.cpp | 7 ++++--- src/spps-agh/calculation_cores/MLTCore.h | 1 + src/spps-agh/data_manager/core_configurationAGH.cpp | 3 +++ src/spps-agh/data_manager/core_configurationAGH.h | 1 + 5 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h index 39c4cfe0..47f44522 100644 --- a/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h +++ b/src/isimpa/data_manager/tree_core/e_core_spps_aghcore_advanced.h @@ -305,6 +305,7 @@ class E_Core_SppsNee_AGH_advanced_MLT : public Element confCore->AppendPropertyBool("selectSeeds", wxTRANSLATE("Use - Seed selection"), true, true); confCore->AppendPropertyDecimal("seedSampleRatio", wxTRANSLATE("Sampling - Seed to sample ratio"), 0.1, false, 2, true, true, 1, 0.01, true); + confCore->AppendPropertyDecimal("mutatedWeight", wxTRANSLATE("Sampling - Mutated particles weight"), 1, false, 1, true, true, 20, 0.5, true); confCore->AppendPropertyDecimal("reflectionPenalty", wxTRANSLATE("Tune - Relfection penalty multiplier"), 0.8, false, 2, true, true, 10, 0.01, true); confCore->AppendPropertyInteger("mutationGeneratorS1", wxTRANSLATE("Tune - Mutation generator S1"), 1024, true, false, true, 0, 1); confCore->AppendPropertyInteger("mutationGeneratorS2", wxTRANSLATE("Tune - Mutation generator S2"), 64, true, false, true, 0, 1); diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index a46df948..a5f287d8 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -17,6 +17,7 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT reflectionPenaltyMultip = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_REFL_PENALTY_MULTIP); seedSampleRatio = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SAMPLE_RATIO); seedSelectionExponent = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_SEED_SELECTION_EXPONENT); + mutatedWeight = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::FPROP_MLT_MUTATED_WEIGHT); selectSeeds = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_SELECT_SEEDS); rejectionSampling = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_DO_REJECTION_SAMPLING); @@ -53,7 +54,7 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) // The ray will be casted as it comes from random distribution, the weight is assigned so that it caries right energy // ((double)InitialSeedsSize/SeedsSize) is essentialy to cancel with weight assigment in CastShadowRay function - inputParticle.weight = 1./ (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize/SeedsSize)); + inputParticle.weight = 1./ (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize/SeedsSize)); RunInitialSeed(inputParticle); if (rejectionSampling) @@ -209,8 +210,8 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu // Treat sample that is not large step if (rejectionSampling){ // Rejection sampling - inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio)) / (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize / SeedsSize)) ; - mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio) + mutatedParticle.largeStep) / (mutation_number * (1. + pLarge) + ((double)InitialSeedsSize / SeedsSize)); + inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio))*mutatedWeight / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)) ; + mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio)*mutatedWeight + mutatedParticle.largeStep) / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)); } else{ // No rejection sampling diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index bff2512b..686e73db 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -24,6 +24,7 @@ class MLTCore : public NextEventEstimationCore float reflectionPenaltyMultip = 0.95; float seedSampleRatio = 0; float seedSelectionExponent = 0.5; + float mutatedWeight = 1.; bool selectSeeds = true; bool rejectionSampling = true; diff --git a/src/spps-agh/data_manager/core_configurationAGH.cpp b/src/spps-agh/data_manager/core_configurationAGH.cpp index 5ed1ce2b..b5cd7436 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.cpp +++ b/src/spps-agh/data_manager/core_configurationAGH.cpp @@ -218,6 +218,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { int mutationGeneratorS1 = advancedNode->GetProperty("mutationGeneratorS1").ToInt(); int mutationGeneratorS2 = advancedNode->GetProperty("mutationGeneratorS2").ToInt(); float mutationProbability = advancedNode->GetProperty("mutationProbability").ToFloat(); + float mutatedWeight = advancedNode->GetProperty("mutatedWeight").ToFloat(); int lowMemoryMode = advancedNode->GetProperty("lowMemoryMode").ToInt(); @@ -229,6 +230,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { std::cout << "Seed selection: " << selectSeeds << std::endl; if (selectSeeds) std::cout << " Seed to sample rato:" << seedSampleRatio << std::endl; if (selectSeeds) std::cout << " Seed selection exponent:" << seedSelectionExponent << std::endl; + std::cout << "Mutated particles weight: " << mutatedWeight << std::endl; std::cout << "Weight Decay Compensation: " << useWeightDecayCompensation << std::endl; std::cout << "Source ray direction mutation: " << mutateSource << std::endl << " Limited functionality, omni source only!!" << std::endl; @@ -268,6 +270,7 @@ void Core_ConfigurationAGH::LoadAdvancedMLT(CXmlNode* simuNode) { SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S1, mutationGeneratorS1); SetConfigInformation(IPROP_MLT_MUTATION_GENERATOR_S2, mutationGeneratorS2); SetConfigInformation(FPROP_MLT_MUTATION_PROBABILITY, mutationProbability); + SetConfigInformation(FPROP_MLT_MUTATED_WEIGHT, mutatedWeight); SetConfigInformation(IPROP_MLT_LOW_MEMORY_MODE, lowMemoryMode); } diff --git a/src/spps-agh/data_manager/core_configurationAGH.h b/src/spps-agh/data_manager/core_configurationAGH.h index b0048c94..65aca256 100644 --- a/src/spps-agh/data_manager/core_configurationAGH.h +++ b/src/spps-agh/data_manager/core_configurationAGH.h @@ -34,6 +34,7 @@ class Core_ConfigurationAGH : public Core_Configuration FPROP_MLT_REFL_PENALTY_MULTIP, FPROP_MLT_SEED_SELECTION_EXPONENT, FPROP_MLT_MUTATION_PROBABILITY, + FPROP_MLT_MUTATED_WEIGHT, }; enum NEW_IPROP From 84b3200a30fe5471e1d9c34d8a3c64cdf7d9b37a Mon Sep 17 00:00:00 2001 From: wbinek Date: Sat, 8 Jan 2022 11:02:57 +0100 Subject: [PATCH 52/56] Increment spps-agh to 0.5.2 --- src/spps-agh/sppsNeeAGHTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index 547cd4b6..b92d757d 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -24,7 +24,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.1beta (23.11.2021). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.2 (08.11.2022). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION From 59501090b58745f11f3030f26dfb3531a50710f5 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 18 Aug 2022 21:39:43 +0200 Subject: [PATCH 53/56] Fix for loading materials from CATT acoustics file --- src/isimpa/data_manager/projet.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/isimpa/data_manager/projet.cpp b/src/isimpa/data_manager/projet.cpp index 88f0938e..b8479ce6 100644 --- a/src/isimpa/data_manager/projet.cpp +++ b/src/isimpa/data_manager/projet.cpp @@ -1054,8 +1054,8 @@ void ProjectManager::ImportMaterialCatt(wxString fromFile) //Attention utilisations d'expressions régulières ( mal de tête en perspective ), mot clé: regex //Fonctionne quelque soit le nombre de bande de frequence - wxRegEx generalMatHeader(wxT("(.*) = <((( ?[+-]?([0-9]*[.])?[0-9]+) ?)*)>( L )?(<((( ?[+-]?([0-9]*[.])?[0-9]+) ?)*)>)?( {([0-9]+) ([0-9]+) ([0-9]+)})?")); - + //wxRegEx generalMatHeader(wxT("(.*) = <((( ?[+-]?([0-9]*[.])?[0-9]+) ?)*)>( L )?(<((( ?[+-]?([0-9]*[.])?[0-9]+) ?)*)>)?( {([0-9]+) ([0-9]+) ([0-9]+)})?")); + wxRegEx generalMatHeader(wxT("(?:ABS|abs|Abs) (.*?)[[:space:]]*=?[[:space:]]*<((( ?[+-]?([0-9]*[.])?[0-9]+) ?)*)>[[:space:]]*(L1?)?[[:space:]]*(<(((.?[+-]?([0-9]*[.])?[0-9]+) ?)*)>)?[[:space:]]*({.?([0-9]+) ([0-9]+) ([0-9]+)})?")); wxRegEx matInfos(wxT("(.*):(.*)")); wxString ligne=lecteur.GetLine(); @@ -1081,7 +1081,7 @@ void ProjectManager::ImportMaterialCatt(wxString fromFile) } wxString absName= generalMatHeader.GetMatch(ligne,1); wxString chaineAbsorption=generalMatHeader.GetMatch(ligne,2); - wxString chaineDiffusion=generalMatHeader.GetMatch(ligne,6); + wxString chaineDiffusion=generalMatHeader.GetMatch(ligne,8); wxString red = generalMatHeader.GetMatch(ligne, 13); wxString green = generalMatHeader.GetMatch(ligne, 14); wxString blue = generalMatHeader.GetMatch(ligne, 15); From 11821fbc3e20ebd50f1660dcc900dd4cffd8ae72 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 20 Apr 2023 02:29:02 +0200 Subject: [PATCH 54/56] Fixes for energy conservation issues, improved sample weighting --- src/spps-agh/calculation_cores/MLTCore.cpp | 71 +++++++++++++++++----- src/spps-agh/calculation_cores/MLTCore.h | 5 +- 2 files changed, 60 insertions(+), 16 deletions(-) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index a5f287d8..50652016 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -36,7 +36,7 @@ MLTCore::MLTCore(t_Mesh& _sceneMesh, t_TetraMesh& _sceneTetraMesh, CONF_CALCULAT lowMemoryMode = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_MLT_LOW_MEMORY_MODE); InitialSeedsSize = *configurationTool->FastGetConfigValue(Core_ConfigurationAGH::IPROP_QUANT_PARTICLE_CALCULATION); - SeedsSize = selectSeeds ? (int)(InitialSeedsSize * seedSampleRatio) : InitialSeedsSize; + TargetSeedsSize = selectSeeds ? (int)(InitialSeedsSize * seedSampleRatio) : InitialSeedsSize; if (mutation_number == 0) { lowMemoryMode = true; @@ -54,7 +54,19 @@ bool MLTCore::Run(CONF_PARTICULE_AGH configurationP) // The ray will be casted as it comes from random distribution, the weight is assigned so that it caries right energy // ((double)InitialSeedsSize/SeedsSize) is essentialy to cancel with weight assigment in CastShadowRay function - inputParticle.weight = 1./ (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize/SeedsSize)); + //float sampleWeightFactor = (double)InitialSeedsSize/(InitialSeedsSize + SeedsSize * mutation_number * (1. * mutatedWeight + pLarge)) * (double)SeedsSize / InitialSeedsSize; + //float sampleWeightFactor = 1. / ((double)InitialSeedsSize / SeedsSize + mutation_number * (1. * mutatedWeight + pLarge)); + + //Wariant 2 + //float sampleWeightFactor = 1. / (mutation_number * (double)SeedsSize / InitialSeedsSize * (1. * mutatedWeight + pLarge) + 1.); + + //Wariant 3 + //float sampleWeightFactor = 1./3.; + + //inputParticle.weight = 1. / (mutation_number * (1. * mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)); + + //Wariant 4 + inputParticle.weight = (double)InitialSeedsSize / ((double)TargetSeedsSize * mutation_number * mutatedWeight + InitialSeedsSize); RunInitialSeed(inputParticle); if (rejectionSampling) @@ -137,12 +149,12 @@ void MLTCore::ChooseSeeds() { // Release memory of the seed on the go InitialSeeds.pop_back(); } - SeedsSize = Seeds.size(); - float avgRayRepeat = (float)SeedsSize / individualRays; + TrueSeedsSize = Seeds.size(); + float avgRayRepeat = (float)TrueSeedsSize / individualRays; - std::cout << "Final seed count: " << SeedsSize << std::endl; - std::cout << "Individual rays: " << individualRays << std::endl; //both values are slightly wrong - std::cout << "Average ray repeat: " << avgRayRepeat << std::endl; //both values are slightly wrong + std::cout << "Final seed count: " << TrueSeedsSize << std::endl; + std::cout << "Individual rays: " << individualRays << std::endl; + std::cout << "Average ray repeat: " << avgRayRepeat << std::endl; if (avgRayRepeat > 10) { std::cout << "WARNING Average ray repeated " << avgRayRepeat <<" times!!" << std::endl << "Consider lower sample to seed ratio" << std::endl; @@ -178,7 +190,7 @@ void MLTCore::CastShadowRays(CONF_PARTICULE_MLT& input_particle, bool delete_sr { //applicationTools.outputTool->NewParticule(confPart); shadowRay.energie *= input_particle.weight; - shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection, if no seed selection the sizes are equal + //shadowRay.energie *= (double)InitialSeedsSize / SeedsSize; // Compensate energy for seed selection, if no seed selection the sizes are equal Run(shadowRay); //applicationTools.outputTool->SaveParticule(); } @@ -190,6 +202,13 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu CONF_PARTICULE_MLT mutatedParticle; double ratio; int mutation_counter = 0; + float mlt_weight; + if (rejectionSampling) { + mlt_weight = ((double)TargetSeedsSize * mutation_number * mutatedWeight) / ((double)TargetSeedsSize * mutation_number * mutatedWeight + InitialSeedsSize); + } + else { + mlt_weight = 1.; + } // We assure that thatu mutation number is greater than 0 in RunMutation function, further logic changed when seed casting for initial seeds was added, to be removed in future //if (mutation_number == 0) @@ -210,27 +229,51 @@ void MLTCore::DoMutations(CONF_PARTICULE_MLT& inputParticle, int totalParticleNu // Treat sample that is not large step if (rejectionSampling){ // Rejection sampling - inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio))*mutatedWeight / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)) ; - mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio)*mutatedWeight + mutatedParticle.largeStep) / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)); - } + + //inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio))*mutatedWeight / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)) ; + //mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio)*mutatedWeight + mutatedParticle.largeStep) / (mutation_number * (1.*mutatedWeight + pLarge) + ((double)InitialSeedsSize / SeedsSize)); + + + //float sampleWeightFactor = (double)SeedsSize / (InitialSeedsSize + SeedsSize * mutation_number * (1. * mutatedWeight + pLarge)); + //float sampleWeightFactor = 1. / ((double)InitialSeedsSize/SeedsSize + mutation_number * (1. * mutatedWeight + pLarge)); + //inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio)) * mutatedWeight * sampleWeightFactor ; + //mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio)*mutatedWeight + mutatedParticle.largeStep) * sampleWeightFactor; + + + // Wariant z wyprowadzeñ z kartki - balance heuristic - wypr2 + //inputParticle.weight += (1. - ratio) * mutatedWeight / (mutation_number * (double)SeedsSize / InitialSeedsSize * (inputParticle.probability.totalProbability(wDC) / b * mutatedWeight + pLarge) + 1); + //mutatedParticle.weight += (ratio * mutatedWeight + mutatedParticle.largeStep) / (mutation_number * (double)SeedsSize / InitialSeedsSize * (mutatedParticle.probability.totalProbability(wDC) / b * mutatedWeight + pLarge) + 1); + + // Wariant z wyprowadzeñ z kartki z wyrzuconym b/I - wypr3 testowo - œrednia z 3 estymatorów + //inputParticle.weight += (b / inputParticle.probability.totalProbability(wDC) * (1. - ratio)) * mutatedWeight / (mutation_number * (double)SeedsSize / InitialSeedsSize * (1 * mutatedWeight + pLarge) + 1); + //mutatedParticle.weight += ((b / mutatedParticle.probability.totalProbability(wDC) * ratio) * mutatedWeight + mutatedParticle.largeStep) / (mutation_number * (double)SeedsSize / InitialSeedsSize * (1 * mutatedWeight + pLarge) + 1); + + // Wariant z wyprowadzeñ z kartki - balance heuristic - wypr4 + inputParticle.weight += (1. - ratio) / (mutation_number * ((double)TrueSeedsSize / InitialSeedsSize) * (inputParticle.probability.totalProbability(wDC) / b + pLarge)); + mutatedParticle.weight += (ratio + mutatedParticle.largeStep) / (mutation_number * ((double)TrueSeedsSize / InitialSeedsSize) * (mutatedParticle.probability.totalProbability(wDC) / b + pLarge)); + } else{ // No rejection sampling - inputParticle.weight = b / inputParticle.probability.totalProbability(wDC) / (1. + pLarge); - mutatedParticle.weight = (b / mutatedParticle.probability.totalProbability(wDC) + mutatedParticle.largeStep) / (1. + pLarge); + inputParticle.weight = b / (inputParticle.probability.totalProbability(wDC) * ((double)TrueSeedsSize / InitialSeedsSize)); + mutatedParticle.weight = b / (mutatedParticle.probability.totalProbability(wDC) * ((double)TrueSeedsSize / InitialSeedsSize)); } //// Make particle seletion, update and SR casting mutatedParticle.largeStep = 0; if(ratio > GetRandValue()) { + inputParticle.weight *= mlt_weight; if (rejectionSampling) CastShadowRays(inputParticle); inputParticle = mutatedParticle; }else - { + { + mutatedParticle.weight *= mlt_weight; if (rejectionSampling) CastShadowRays(mutatedParticle); } mutation_counter++; } + + inputParticle.weight *= mlt_weight; CastShadowRays(inputParticle); } diff --git a/src/spps-agh/calculation_cores/MLTCore.h b/src/spps-agh/calculation_cores/MLTCore.h index 686e73db..3a226c72 100644 --- a/src/spps-agh/calculation_cores/MLTCore.h +++ b/src/spps-agh/calculation_cores/MLTCore.h @@ -46,8 +46,9 @@ class MLTCore : public NextEventEstimationCore std::list InitialSeeds; std::list Seeds; - int InitialSeedsSize=0; - int SeedsSize=0; + int InitialSeedsSize = 0; + int TargetSeedsSize = 0; + int TrueSeedsSize = 0; void Mutate(const CONF_PARTICULE_MLT& input_particle, CONF_PARTICULE_MLT& mutetedParticle) const; bool MoveToNextReflection(CONF_PARTICULE_MLT& configurationP, double rndDir1, double rndDir2, double rndTravelProbability, double rndMatAbsorbtion, float deltaT, float distanceToTravel); From 9872a2164e5a5df510eae4c3e1448f07a849f2de Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 20 Apr 2023 02:29:19 +0200 Subject: [PATCH 55/56] Increment SPPS to 0.5.3 --- src/spps-agh/sppsNeeAGHTypes.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/spps-agh/sppsNeeAGHTypes.h b/src/spps-agh/sppsNeeAGHTypes.h index b92d757d..e1214ecd 100644 --- a/src/spps-agh/sppsNeeAGHTypes.h +++ b/src/spps-agh/sppsNeeAGHTypes.h @@ -24,7 +24,7 @@ -#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.2 (08.11.2022). Based on SPPS Nantes 2.2.1 version ocober 2020" +#define SPPSNEE_AGH_VERSION "SPPS-AGH 0.5.3 (19.04.2023). Based on SPPS Nantes 2.2.1 version ocober 2020" #define __USE_BOOST_RANDOM_GENERATOR__ #define UTILISER_MAILLAGE_OPTIMISATION @@ -64,7 +64,7 @@ struct CONF_PARTICULE_MLT : public CONF_PARTICULE_AGH int largeStep = 0; bool isInitialSeed = false; - class + class PartProbability { public: double currentProbability = 1; //probability of generating ray path up to some point From c991e13ef24eabd68ff16eb31db83b9660d97e82 Mon Sep 17 00:00:00 2001 From: wbinek Date: Thu, 20 Apr 2023 14:00:59 +0200 Subject: [PATCH 56/56] Regression, MLT calculation with seed selection disabled --- src/spps-agh/calculation_cores/MLTCore.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spps-agh/calculation_cores/MLTCore.cpp b/src/spps-agh/calculation_cores/MLTCore.cpp index 50652016..b861e3bb 100644 --- a/src/spps-agh/calculation_cores/MLTCore.cpp +++ b/src/spps-agh/calculation_cores/MLTCore.cpp @@ -111,6 +111,7 @@ void MLTCore::ChooseSeeds() { // If seed chosing disabled if (!selectSeeds) { Seeds = std::move(InitialSeeds); + TrueSeedsSize = Seeds.size(); std::cout << "Seed selection disabled, copied all initial seeds" << std::endl; return; }