Skip to content

Commit

Permalink
Add functions for calculating coefficients
Browse files Browse the repository at this point in the history
  • Loading branch information
Dylan Harries committed Jan 19, 2017
1 parent 7c8758f commit e10455c
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 23 deletions.
11 changes: 8 additions & 3 deletions meta/FlexibleSUSY.m
Expand Up @@ -1459,23 +1459,28 @@ corresponding tadpole is real or imaginary (only in models with CP

WriteSemiAnalyticSolutionsClass[semiAnalyticBCs_List, semiAnalyticSolns_List, files_List] :=
Module[{semiAnalyticSolutionsDefs = "", boundaryValueStructDefs = "",
coefficientGetters = "", basisLambdas = "", applyBoundaryConditions = "",
coefficientGetters = "", createBasisEvaluators = "", applyBoundaryConditions = "",
datasets, numberOfTrialPoints, initializeTrialBoundaryValues = "",
createLinearSystemSolvers = "", calculateCoefficients = "",
calculateCoefficientsPrototypes = "", calculateCoefficientsFunctions = ""},
semiAnalyticSolutionsDefs = SemiAnalytic`CreateSemiAnalyticSolutionsDefinitions[semiAnalyticSolns];
boundaryValueStructDefs = SemiAnalytic`CreateLocalBoundaryValuesDefinitions[semiAnalyticSolns];
coefficientGetters = SemiAnalytic`CreateSemiAnalyticCoefficientGetters[semiAnalyticSolns];
basisLambdas = SemiAnalytic`CreateBasisEvaluators[semiAnalyticSolns];
createBasisEvaluators = SemiAnalytic`CreateBasisEvaluators[semiAnalyticSolns];
datasets = SemiAnalytic`ConstructTrialDatasets[semiAnalyticSolns];
createLinearSystemSolvers = SemiAnalytic`CreateLinearSystemSolvers[datasets, semiAnalyticSolns];
{numberOfTrialPoints, initializeTrialBoundaryValues} = SemiAnalytic`InitializeTrialInputValues[datasets];
applyBoundaryConditions = SemiAnalytic`ApplySemiAnalyticBoundaryConditions[semiAnalyticBCs, semiAnalyticSolns];
calculateCoefficients = SemiAnalytic`CalculateCoefficients[datasets];
{calculateCoefficientsPrototypes, calculateCoefficientsFunctions} = SemiAnalytic`CreateCoefficientsCalculations[semiAnalyticSolns];
WriteOut`ReplaceInFiles[files, { "@semiAnalyticSolutionsDefs@" -> IndentText[WrapLines[semiAnalyticSolutionsDefs]],
"@boundaryValueStructDefs@" -> IndentText[IndentText[WrapLines[boundaryValueStructDefs]]],
"@coefficientGetters@" -> IndentText[WrapLines[coefficientGetters]],
"@numberOfTrialPoints@" -> ToString[numberOfTrialPoints],
"@initializeTrialBoundaryValues@" -> IndentText[WrapLines[initializeTrialBoundaryValues]],
"@basisLambdas@" -> IndentText[WrapLines[basisLambdas]],
"@createBasisEvaluators@" -> IndentText[WrapLines[createBasisEvaluators]],
"@createLinearSystemSolvers@" -> IndentText[WrapLines[createLinearSystemSolvers]],
"@calculateCoefficients@" -> IndentText[calculateCoefficients],
"@applyBoundaryConditions@" -> IndentText[WrapLines[applyBoundaryConditions]],
"@calculateCoefficientsPrototypes@" -> IndentText[calculateCoefficientsPrototypes],
"@calculateCoefficientsFunctions@" -> calculateCoefficientsFunctions,
Expand Down
47 changes: 43 additions & 4 deletions meta/SemiAnalytic.m
Expand Up @@ -39,6 +39,8 @@
{integer id, {pars}, {input values}}";
InitializeTrialInputValues::usage="";
CreateBasisEvaluators::usage="";
CreateLinearSystemSolvers::usage="";
CalculateCoefficients::usage="";
CreateSemiAnalyticCoefficientsCalculation::usage="";
CreateCoefficientsCalculations::usage="";

Expand Down Expand Up @@ -596,10 +598,6 @@
Return[result];
];

CreateSemiAnalyticCoefficientsCalculation[solutions_List] :=
Module[{trialValues},
];

ApplySettingLocally[{parameter_, value_}, modelPrefix_String] :=
Which[Parameters`IsModelParameter[parameter],
Parameters`SetParameter[parameter, value, modelPrefix],
Expand Down Expand Up @@ -728,6 +726,47 @@
Utils`StringJoinWithSeparator[evaluators, "\n"] <> "\n"
];

CreateLinearSystemSolver[dataset_, solutions_List] :=
Module[{idx, inputs, pars, parTypes, elementType, basisLengths, basisSize,
solverName, evaluatorName, result = ""},
idx = ToString[dataset[[1]]];
pars = dataset[[2]];
parTypes = CConversion`GetScalarElementType[Parameters`GetType[#]]& /@ dataset[[2]];
If[MemberQ[parTypes, CConversion`ScalarType[CConversion`complexScalarCType]],
elementType = CConversion`CreateCType[CConversion`ScalarType[CConversion`complexScalarCType]],
elementType = CConversion`CreateCType[CConversion`ScalarType[CConversion`realScalarCType]]
];
basisLengths = Length[GetBasis[#]]& /@ Select[solutions, MemberQ[pars, GetName[#]]&];
If[Length[DeleteDuplicates[basisLengths]] =!= 1,
Print["Error: invalid collection of parameters specified:", pars];
Quit[1];
];
basisSize = ToString[First[basisLengths]];
solverName = "solver_" <> idx;
evaluatorName = "basis_" <> idx;
"auto " <> solverName <> " = create_solver<" <> elementType
<> "," <> basisSize <> ">(datasets[" <> idx <> "], " <> evaluatorName <> ");\n"
];

CreateLinearSystemSolvers[datasets_List, solutions_List] :=
Module[{result = ""},
(result = result <> CreateLinearSystemSolver[#, solutions])& /@ datasets;
Return[result];
];


CalculateCoefficients[datasets_List] :=
Module[{i, result = ""},
For[i = 1, i <= Length[datasets], i++,
(result = result <> "calculate_"
<> CConversion`ToValidCSymbolString[#]
<> "_coefficients(solver_" <> ToString[datasets[[i,1]]]
<> ", datasets[" <> ToString[datasets[[i,1]]]
<> "]);\n")& /@ datasets[[i,2]];
];
Return[result];
];

EvaluateSemiAnalyticSolution[solution_] :=
Module[{parameter, basisRules, coeffs},
parameter = GetName[solution];
Expand Down
8 changes: 4 additions & 4 deletions templates/semi_analytic_solutions.cpp.in
Expand Up @@ -72,7 +72,7 @@ std::map<int,@ModelName@_semi_analytic_solutions::Data_vector_t> @ModelName@_sem
std::map<int,Data_vector_t> datasets;
for (const auto& point: trial_data) {
for (auto idx: point.basis_sets) {
datasets[idx].push_back(&point.model);
datasets[idx].push_back(&point);
}
}
return datasets;
Expand All @@ -83,11 +83,11 @@ void @ModelName@_semi_analytic_solutions::calculate_coefficients(const @ModelNam
@ModelName@_soft_parameters model(model_);
model.run_to(input_scale);

@basisLambdas@

calculate_trial_data(model);
std::map<int,Data_vector_t> datasets(create_datasets());

// for each parameter, compute coefficients using appropriate dataset
@createBasisEvaluators@
@createLinearSystemSolvers@
@calculateCoefficients@
}

Expand Down
26 changes: 14 additions & 12 deletions templates/semi_analytic_solutions.hpp.in
Expand Up @@ -60,8 +60,6 @@ public:
void calculate_coefficients(const @ModelName@_soft_parameters&);

private:
using Data_vector_t = std::vector<@ModelName@_soft_parameters*>

struct Boundary_values {
@boundaryValueStructDefs@
};
Expand All @@ -72,6 +70,8 @@ private:
std::vector<int> basis_sets{};
};

using Data_vector_t = std::vector<Model_data*>

double input_scale{0.}; ///< scale at which boundary conditions hold
double output_scale{0.}; ///< scale at which coefficients are calculated
std::array<Model_data,@numberOfTrialPoints@> trial_data{};
Expand All @@ -84,22 +84,24 @@ private:
void set_to_boundary_values(@ModelName@_soft_parameters&, const Boundary_values&) const;
std::map<int,Data_vector_t> create_datasets() const;

template <class Scalar, int Rows, int Cols, class Basis>
Eigen::Matrix<Scalar,Rows,Cols> compute_basis_points(const std::vector<Boundary_values>&, const Basis&) const;

template <class Scalar, int BasisSize, class BasisEvaluator>
Eigen::JacobiSVD<Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> > create_solver(
const Data_vector_t&, const BasisEvaluator&) const;

@calculateCoefficientsPrototypes@
};

template <class Scalar, int Rows, int Cols, class Basis>
Eigen::Matrix<Scalar,Rows,Cols> @ModelName@_semi_analytic_solutions::compute_basis_points(
const std::vector<Boundary_values>& points, const Basis& basis) const
template <class Scalar, int BasisSize, class BasisEvaluator>
Eigen::JacobiSVD<Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> > create_solver(
const Data_vector_t& data, const BasisEvaluator& ev) const
{
Eigen::Matrix<Scalar,Rows,Cols> result;
for (int i = 0; i < Rows; ++i) {
result.row(i) = basis(points.at(i));
const std::size_t n = data.size();
Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> lhs(n, BasisSize);
for (std::size_t i = 0; i < n; ++i) {
lhs.row(i) = ev(data[i]->boundary_values);
}
return result;
return Eigen::JacobiSVD<Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> >(
lhs, Eigen::ComputeThinU | Eigen::ComputeThinV);
}

} // namespace flexiblesusy
Expand Down

0 comments on commit e10455c

Please sign in to comment.