CDPP / AMDA_Kernel

Commit c34a11426f88e0c9a289e96663ac807a43cbdf47

Authored by Benjamin Renard
2 parents 83cbfd9f e8165c26
Exists in master and in 37 other branches AllLegends, FER_10400, FER_10412, FER_11000, FER_11100, FER_11242, FER_11277, FER_11500, FER_11537, FER_11588, FER_11651, FER_19836, FER_5344, FER_6031, FER_7318, FER_7618, FER_7857, FER_8455, FER_9629, FER_CR41, FER_CR42, FER_CR43, FER_TFCat, Histo1D2D, JUNO_JEDI, MAZ_10936, MAZ_6018, MAZ_9016_New, MAZ_FIX_9016, SpeasyGet, amdadev, cu2022, das2, density_maps, develop, evol_9042, ubuntu

Merge branch 'develop' into MAZ_6018

Showing 8 changed files with 392 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/ExternLib/Tetrahedron/common/tetrahedron.cc 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,180 @@
  1 +#include<vector>
  2 +#include<iostream>
  3 +#include<cmath>
  4 +#include<algorithm>
  5 +
  6 +using Vector = std::vector<double>;
  7 +using Matrix = std::vector<Vector>;
  8 +
  9 +#ifdef DEBUG
  10 +void print_matrix(Matrix);
  11 +#endif
  12 +
  13 +Matrix concatenate_vectors(Vector p1, Vector p2, Vector p3, Vector p4) {
  14 + Matrix mat;
  15 + mat.push_back(p1);
  16 + mat.push_back(p2);
  17 + mat.push_back(p3);
  18 + mat.push_back(p4);
  19 + return mat;
  20 +}
  21 +
  22 +Matrix transpose(Matrix mat) {
  23 + Matrix t; // new matrix
  24 + int n = mat.size();
  25 + int m = mat[0].size();
  26 + for(int i=0;i<m;i++) {
  27 + // new row of size n
  28 + Vector r;
  29 + for(int j=0;j<n;j++) {
  30 + r.push_back(mat[j][i]);
  31 + }
  32 + t.push_back(r);
  33 + }
  34 + return t;
  35 +}
  36 +
  37 +double scalar_product(Vector vect1, Vector vect2) {
  38 + double ans=0;
  39 + for(int i=0;i<vect1.size();i++) {
  40 + ans += vect1[i] * vect2[i];
  41 + }
  42 + return ans;
  43 +}
  44 +
  45 +Matrix matrix_multiply(Matrix mat1, Matrix mat2) {
  46 + Matrix mat2t = transpose(mat2);
  47 + Matrix r;
  48 + int n = mat1.size();
  49 + int m = mat1[0].size();
  50 + for(int i=0;i<n;i++) {
  51 + // new row
  52 + Vector row;
  53 + for(int j=0;j<n;j++) row.push_back(scalar_product(mat1[i], mat2t[j]));
  54 + r.push_back(row);
  55 + }
  56 + return r;
  57 +
  58 +}
  59 +
  60 +#ifdef DEBUG
  61 +void print_vector(Vector vect) {
  62 + for(int i=0;i<vect.size();i++) std::cout << vect[i] << " ";
  63 + std::cout << std::endl;
  64 +}
  65 +
  66 +void print_matrix(Matrix mat) {
  67 + for(int i=0;i<mat.size();i++) {
  68 + print_vector(mat[i]);
  69 + }
  70 +}
  71 +#endif
  72 +
  73 +Vector vector_substract(Vector vect1, Vector vect2) {
  74 + Vector ans;
  75 + for(int i=0;i<vect1.size();i++) {
  76 + ans.push_back(vect1[i] - vect2[i]);
  77 + }
  78 + return ans;
  79 +}
  80 +
  81 +Vector vector_projection(Vector v, Vector u) {
  82 + Vector ans;
  83 + double a = scalar_product(u,v);
  84 + double b = scalar_product(u,u);
  85 + for(int i=0;i<v.size();i++) {
  86 + ans.push_back((a/b) * u[i]);
  87 + }
  88 + return ans;
  89 +}
  90 +
  91 +Vector vector_normalize(Vector vect) {
  92 + Vector ans;
  93 + double norm = scalar_product(vect, vect);
  94 + for(int i=0;i<vect.size();i++) ans.push_back(vect[i] / sqrt(norm));
  95 + return ans;
  96 +}
  97 +
  98 +Matrix gram_schmidt(Matrix mat) {
  99 + Matrix matt = transpose(mat);
  100 + Matrix basis;
  101 + Matrix Us;
  102 + for(int i=0;i<mat.size();i++) {
  103 + // new basis vector
  104 + Vector nb = matt[i];
  105 + for(int j=0;j<=i-1;j++) {
  106 + nb = vector_substract(nb, vector_projection(matt[i], Us[j]));
  107 + Vector vv = vector_projection(matt[i], Us[j]);
  108 + }
  109 + Us.push_back(nb);
  110 + nb = vector_normalize(nb);
  111 + basis.push_back(vector_normalize(nb));
  112 + }
  113 + return transpose(basis);
  114 +}
  115 +
  116 +bool is_upper_triangular(Matrix mat) {
  117 + for(int i=0;i<mat.size();i++) {
  118 + for(int j=0;j<i;j++) {
  119 + if(mat[i][j] > 1e-16) return false;
  120 + }
  121 + }
  122 + return true;
  123 +}
  124 +
  125 +Vector matrix_eigenvalues(Matrix mat) {
  126 + Matrix A = mat;
  127 + for(int k=0;k<1000000;k++) {
  128 + Matrix basis = gram_schmidt(A);
  129 + Matrix basisT = transpose(basis);
  130 + Matrix t1 = matrix_multiply(A, basis);
  131 + Matrix t2 = matrix_multiply(basisT, t1);
  132 + if(is_upper_triangular(t2)) break;
  133 +
  134 + A = t2;
  135 + }
  136 + Vector eig;
  137 + for(int i=0;i<mat.size();i++) eig.push_back(A[i][i]);
  138 + return eig;
  139 +}
  140 +
  141 +void compute_tetrahedron(Vector p1, Vector p2, Vector p3, Vector p4, double& elongation, double& planarity) {
  142 + // make matrix by concatenating vectors
  143 + Matrix mat = concatenate_vectors(p1,p2,p3,p4);
  144 +
  145 + Matrix matt = transpose(mat);
  146 +
  147 + // matrix product
  148 + Matrix pp = matrix_multiply(matt, mat);
  149 +
  150 + Vector eig = matrix_eigenvalues(pp);
  151 +
  152 + // sort eigenvalues
  153 + sort(eig.begin(), eig.end());
  154 +
  155 + elongation = 1. - sqrt(eig[0] / eig[1]);
  156 + planarity = 1. - sqrt(eig[1] / eig[2]);
  157 +}
  158 +
  159 +#ifdef DEBUG
  160 +Vector random_vector() {
  161 + Vector v;
  162 + for(int i=0;i<3;i++) v.push_back((double)(rand()%7));
  163 +
  164 + return v;
  165 +}
  166 +
  167 +int main(int argc, char ** argv) {
  168 + // create some random data and compute the elongation
  169 + for(int i=0;i<1000000;i++) {
  170 + // generate twelve random numbers
  171 + Vector p1 = random_vector();
  172 + Vector p2 = random_vector();
  173 + Vector p3 = random_vector();
  174 + Vector p4 = random_vector();
  175 + std::cout << "Elongation " << elongation(p1,p2,p3,p4) << ", planarity: " << planarity(p1,p2,p3,p4) << std::endl;
  176 + }
  177 + return 0;
  178 +}
  179 +#endif
  180 +
... ...
src/ExternLib/Tetrahedron/common/tetrahedron.hh 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,28 @@
  1 +/**
  2 + * tetrahedron.hh
  3 + *
  4 + * Created on: 12 jul. 2022
  5 + * Author: Alexandre Schulz - BRE
  6 + */
  7 +
  8 +#ifndef TETRAHEDRON_HH_
  9 +#define TETRAHEDRON_HH_
  10 +
  11 +#include <vector>
  12 +
  13 +using Vector = std::vector<double>;
  14 +using Matrix = std::vector<Vector>;
  15 +
  16 +void compute_tetrahedron(Vector p1, Vector p2, Vector p3, Vector p4, double& elongation, double& planarity);
  17 +
  18 +template<typename Type>
  19 +void compute_tetrahedron_(std::vector<Type> p1, std::vector<Type> p2, std::vector<Type> p3, std::vector<Type> p4, double& elongation, double& planarity) {
  20 + std::vector<double> d_p1(p1.begin(), p1.end());
  21 + std::vector<double> d_p2(p2.begin(), p2.end());
  22 + std::vector<double> d_p3(p3.begin(), p3.end());
  23 + std::vector<double> d_p4(p4.begin(), p4.end());
  24 + compute_tetrahedron(d_p1, d_p2, d_p3, d_p4, elongation, planarity);
  25 +}
  26 +
  27 +#endif /* TETRAHEDRON_HH_ */
  28 +
... ...
src/ExternLib/Tetrahedron/elongation/CMakeLists.txt 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,25 @@
  1 +
  2 +PROJECT(elongation)
  3 +
  4 +set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH}/elongation)
  5 +
  6 +include_directories(
  7 + ../common
  8 +)
  9 +
  10 +#Library configuration
  11 +file(
  12 + GLOB_RECURSE
  13 + source_files
  14 + ./*
  15 + ../common/*
  16 +)
  17 +
  18 +ADD_LIBRARY( elongation SHARED ${source_files} )
  19 +
  20 +target_link_libraries(
  21 + elongation
  22 +)
  23 +
  24 +FILE( GLOB_RECURSE PROJ_HEADERS *.hh )
  25 +INSTALL(FILES ${PROJ_HEADERS} DESTINATION ${LIBRARY_OUTPUT_PATH} PERMISSIONS OWNER_READ GROUP_READ WORLD_READ)
... ...
src/ExternLib/Tetrahedron/elongation/elongation.cc 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,43 @@
  1 +#include<vector>
  2 +#include<iostream>
  3 +#include<cmath>
  4 +#include<algorithm>
  5 +
  6 +#include "elongation.hh"
  7 +#include "tetrahedron.hh"
  8 +
  9 +double elongation(const std::vector<float>& p1, const std::vector<float>& p2, const std::vector<float>& p3, const std::vector<float>& p4) {
  10 + double elongation = 0.;
  11 + double planarity = 0.;
  12 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  13 + return elongation;
  14 +}
  15 +
  16 +double elongation(const std::vector<double>& p1, const std::vector<double>& p2, const std::vector<double>& p3, const std::vector<double>& p4) {
  17 + double elongation = 0.;
  18 + double planarity = 0.;
  19 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  20 + return elongation;
  21 +}
  22 +
  23 +double elongation(const std::vector<long double>& p1, const std::vector<long double>& p2, const std::vector<long double>& p3, const std::vector<long double>& p4) {
  24 + double elongation = 0.;
  25 + double planarity = 0.;
  26 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  27 + return elongation;
  28 +}
  29 +
  30 +double elongation(const std::vector<int>& p1, const std::vector<int>& p2, const std::vector<int>& p3, const std::vector<int>& p4) {
  31 + double elongation = 0.;
  32 + double planarity = 0.;
  33 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  34 + return elongation;
  35 +}
  36 +
  37 +double elongation(const std::vector<short>& p1, const std::vector<short>& p2, const std::vector<short>& p3, const std::vector<short>& p4) {
  38 + double elongation = 0.;
  39 + double planarity = 0.;
  40 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  41 + return elongation;
  42 +}
  43 +
... ...
src/ExternLib/Tetrahedron/elongation/elongation.hh 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,24 @@
  1 +/**
  2 + * elongation.hh
  3 + *
  4 + * Created on: 12 jul. 2022
  5 + * Author: Alexandre Schulz - BRE
  6 + */
  7 +
  8 +#ifndef ELONGATION_HH_
  9 +#define ELONGATION_HH_
  10 +
  11 +#include <vector>
  12 +
  13 +double elongation(const std::vector<float>& p1, const std::vector<float>& p2, const std::vector<float>& p3, const std::vector<float>& p4);
  14 +
  15 +double elongation(const std::vector<double>& p1, const std::vector<double>& p2, const std::vector<double>& p3, const std::vector<double>& p4);
  16 +
  17 +double elongation(const std::vector<long double>& p1, const std::vector<long double>& p2, const std::vector<long double>& p3, const std::vector<long double>& p4);
  18 +
  19 +double elongation(const std::vector<int>& p1, const std::vector<int>& p2, const std::vector<int>& p3, const std::vector<int>& p4);
  20 +
  21 +double elongation(const std::vector<short>& p1, const std::vector<short>& p2, const std::vector<short>& p3, const std::vector<short>& p4);
  22 +
  23 +#endif /* ELONGATION_HH_ */
  24 +
... ...
src/ExternLib/Tetrahedron/planarity/CMakeLists.txt 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,25 @@
  1 +
  2 +PROJECT(planarity)
  3 +
  4 +set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH}/planarity)
  5 +
  6 +include_directories(
  7 + ../common/
  8 +)
  9 +
  10 +#Library configuration
  11 +file(
  12 + GLOB_RECURSE
  13 + source_files
  14 + ./*
  15 + ../common/*
  16 +)
  17 +
  18 +ADD_LIBRARY( planarity SHARED ${source_files} )
  19 +
  20 +target_link_libraries(
  21 + planarity
  22 +)
  23 +
  24 +FILE( GLOB_RECURSE PROJ_HEADERS *.hh )
  25 +INSTALL(FILES ${PROJ_HEADERS} DESTINATION ${LIBRARY_OUTPUT_PATH} PERMISSIONS OWNER_READ GROUP_READ WORLD_READ)
... ...
src/ExternLib/Tetrahedron/planarity/planarity.cc 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,43 @@
  1 +#include<vector>
  2 +#include<iostream>
  3 +#include<cmath>
  4 +#include<algorithm>
  5 +
  6 +#include "planarity.hh"
  7 +#include "tetrahedron.hh"
  8 +
  9 +double planarity(const std::vector<float>& p1, const std::vector<float>& p2, const std::vector<float>& p3, const std::vector<float>& p4) {
  10 + double elongation = 0.;
  11 + double planarity = 0.;
  12 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  13 + return planarity;
  14 +}
  15 +
  16 +double planarity(const std::vector<double>& p1, const std::vector<double>& p2, const std::vector<double>& p3, const std::vector<double>& p4) {
  17 + double elongation = 0.;
  18 + double planarity = 0.;
  19 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  20 + return planarity;
  21 +}
  22 +
  23 +double planarity(const std::vector<long double>& p1, const std::vector<long double>& p2, const std::vector<long double>& p3, const std::vector<long double>& p4) {
  24 + double elongation = 0.;
  25 + double planarity = 0.;
  26 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  27 + return planarity;
  28 +}
  29 +
  30 +double planarity(const std::vector<int>& p1, const std::vector<int>& p2, const std::vector<int>& p3, const std::vector<int>& p4) {
  31 + double elongation = 0.;
  32 + double planarity = 0.;
  33 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  34 + return planarity;
  35 +}
  36 +
  37 +double planarity(const std::vector<short>& p1, const std::vector<short>& p2, const std::vector<short>& p3, const std::vector<short>& p4) {
  38 + double elongation = 0.;
  39 + double planarity = 0.;
  40 + compute_tetrahedron_(p1, p2, p3, p4, elongation, planarity);
  41 + return planarity;
  42 +}
  43 +
... ...
src/ExternLib/Tetrahedron/planarity/planarity.hh 0 → 100644
  Diff comments   View file @c34a114
... ... @@ -0,0 +1,24 @@
  1 +/**
  2 + * planarity.hh
  3 + *
  4 + * Created on: 12 jul. 2022
  5 + * Author: Alexandre Schulz - BRE
  6 + */
  7 +
  8 +#ifndef PLANARITY_HH_
  9 +#define PLANARITY_HH_
  10 +
  11 +#include <vector>
  12 +
  13 +double planarity(const std::vector<float>& p1, const std::vector<float>& p2, const std::vector<float>& p3, const std::vector<float>& p4);
  14 +
  15 +double planarity(const std::vector<double>& p1, const std::vector<double>& p2, const std::vector<double>& p3, const std::vector<double>& p4);
  16 +
  17 +double planarity(const std::vector<long double>& p1, const std::vector<long double>& p2, const std::vector<long double>& p3, const std::vector<long double>& p4);
  18 +
  19 +double planarity(const std::vector<int>& p1, const std::vector<int>& p2, const std::vector<int>& p3, const std::vector<int>& p4);
  20 +
  21 +double planarity(const std::vector<short>& p1, const std::vector<short>& p2, const std::vector<short>& p3, const std::vector<short>& p4);
  22 +
  23 +#endif /* PLANARITY_HH_ */
  24 +
... ...