Skip to content

Commit b31a79f

Browse files
author
Kun Chen
committed
this commit introduce a bug. need to find it
1 parent 60e305d commit b31a79f

10 files changed

+231
-255
lines changed

‎in00

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
1-
20 1.0 1.0 3.0 6 0 -15 1
1+
20 1.0 1.0 3.0 101 0 -15 1
22

33
# Beta, rs, Mass2, MaxExtMom(*kF), TotalStep(*1e6), Observable, Seed, PID

‎inlist

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
20 1.0 1.0 3.0 101 0 16
1+
20 1.0 1.0 3.0 101 0 8
22

33
10 1.0 1.0 10.0 101 1 8
44
# Beta, rs, Mass2, MaxExtMom(*kF), TotalStep(*1e6), Observable, Duplication

‎polar_freq.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@
4747
Num=0
4848
data0=None
4949
for f in files:
50-
if re.match("group"+str(order-1)+"_pid[0-9]+.dat", f):
50+
if re.match("group"+str(order)+"_pid[0-9]+.dat", f):
5151
print f
5252
Num+=1
5353
d=np.loadtxt(folder+f)

‎src/global.h

+5-4
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
#ifndef FeynCalc_global_h
22
#define FeynCalc_global_h
33

4+
#include "utility/vector.h"
45
#include <math.h>
5-
// #include "utility/vector.h"
6-
#include <array>
6+
#include <vector>
77

88
// turn off all assert
99
const bool DEBUGMODE = true;
@@ -30,6 +30,7 @@ struct parameter {
3030
int SaveFileTimer; // how many secondes between saving to file
3131
int MessageTimer; // how many secondes between two checking for message
3232
int ReweightTimer; // how many secondes between two reweighting
33+
std::vector<int> GroupID;
3334
};
3435

3536
/////////// Global Constants ////////////////////
@@ -79,8 +80,8 @@ const int INL = 0, INR = 1, OUTL = 3, OUTR = 4;
7980
#define FLIP(x) (1 - x)
8081
//////////////////////////////////////////////////////
8182

82-
// typedef Vec<double, D> mom;
83+
typedef Vec<double, D> momentum;
8384

84-
typedef std::array<double, D> momentum;
85+
// typedef std::array<double, D> momentum;
8586

8687
#endif

‎src/main.cpp

+6
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,12 @@ void InitPara() {
7373
Para.PrinterTimer = 5;
7474
Para.SaveFileTimer = 30;
7575
Para.ReweightTimer = 30;
76+
77+
Para.GroupID = {
78+
1,
79+
2,
80+
3,
81+
};
7682
}
7783

7884
void MonteCarlo(markov &Markov) {

‎src/markov.cpp

+29-28
Original file line numberDiff line numberDiff line change
@@ -82,10 +82,8 @@ void markov::Initialization(string FilePrefix) {
8282

8383
///==== initialize observable =======================//
8484
for (auto &g : Weight.Groups) {
85-
polar p = polar();
86-
p.fill(1.0e-10);
87-
Polar.push_back(p);
88-
PolarStatic.push_back(1.0e-10);
85+
Polar[g.ID].fill(1.0e-10);
86+
PolarStatic[g.ID] = 1.0e-10;
8987
}
9088
///=== Do all kinds of test =======================//
9189
Weight.StaticTest();
@@ -104,24 +102,24 @@ void markov::AdjustGroupReWeight() {
104102
void markov::Measure() {
105103
double MCWeight = fabs(Var.CurrGroup->Weight) * Var.CurrGroup->ReWeight;
106104
double WeightFactor = Var.CurrGroup->Weight / MCWeight;
105+
107106
Polar[Var.CurrGroup->ID][Var.CurrExtMomBin] += WeightFactor;
108107
PolarStatic[Var.CurrGroup->ID] += WeightFactor;
109108
};
110109

111110
void markov::SaveToFile() {
112-
for (int i = 0; i < Polar.size(); i++) {
111+
for (auto &id : Para.GroupID) {
113112
ofstream PolarFile;
114-
string FileName =
115-
tfm::format("group%d_pid%d.dat", Weight.Groups[i].ID, Para.PID);
113+
string FileName = tfm::format("group%d_pid%d.dat", id, Para.PID);
116114
PolarFile.open(FileName, ios::out | ios::trunc);
117115
if (PolarFile.is_open()) {
118116
PolarFile << tfm::format(
119-
"#PID: %d, rs:%.3f, Beta: %.3f, Group: %d, Step: %d\n", Para.PID,
120-
Para.Rs, Para.Beta, Weight.Groups[i].ID, Para.Counter);
117+
"#PID: %d, Type:%d, rs:%.3f, Beta: %.3f, Group: %d, Step: %d\n",
118+
Para.PID, Para.ObsType, Para.Rs, Para.Beta, id, Para.Counter);
121119

122-
for (int j = 0; j < Polar[i].size(); j++)
120+
for (int j = 0; j < Polar[id].size(); j++)
123121
PolarFile << tfm::format("%13.6f\t%13.6f\n", Var.ExtMomTable[j][0],
124-
Polar[i][j]);
122+
Polar[id][j]);
125123
PolarFile.close();
126124
} else {
127125
LOG_WARNING("Polarization for PID " << Para.PID << " fails to save!");
@@ -132,12 +130,12 @@ void markov::SaveToFile() {
132130
string FileName = tfm::format("output%d.dat", Para.PID);
133131
StaticPolarFile.open(FileName, ios::out | ios::trunc);
134132
if (StaticPolarFile.is_open()) {
135-
for (int i = 0; i < PolarStatic.size(); i++) {
133+
for (auto &id : Para.GroupID) {
136134
StaticPolarFile << tfm::format(
137-
"PID:%-4d Group:%-4d rs:%-.3f "
135+
"PID:%-4d Type:%-2d Group:%-4d rs:%-.3f "
138136
"Beta:%-.3f Lambda:%-.3f Polar: % 13.6f\n",
139-
Para.PID, Weight.Groups[i].ID, Para.Rs, Para.Beta, Para.Mass2,
140-
PolarStatic[i]);
137+
Para.PID, Para.ObsType, id, Para.Rs, Para.Beta, Para.Mass2,
138+
PolarStatic[id]);
141139
}
142140
StaticPolarFile.close();
143141
} else {
@@ -233,9 +231,6 @@ void markov::ChangeTau() {
233231
Var.Tau[TauIndex] = NewTau;
234232
Var.Tau[TauIndex + 1] = NewTau;
235233

236-
// if (Var.CurrGroup->ID == 1 && Para.Counter > 100000)
237-
// cout << Var.CurrGroup->ID << endl;
238-
239234
Weight.ChangeTau(*Var.CurrGroup, TauIndex);
240235
double NewWeight = Weight.GetNewWeight(*Var.CurrGroup);
241236
double R = Prop * fabs(NewWeight) / fabs(Var.CurrGroup->Weight);
@@ -252,22 +247,26 @@ void markov::ChangeTau() {
252247
};
253248

254249
void markov::ChangeMomentum() {
255-
static momentum CurrMom;
256250
int LoopIndex = Random.irn(0, Var.CurrGroup->LoopNum - 1);
257251
Proposed[CHANGE_MOM][Var.CurrGroup->ID]++;
258252

259-
COPYFROMTO(Var.LoopMom[LoopIndex], CurrMom);
260-
261253
double Prop;
262254
int NewExtMomBin;
255+
static momentum CurrMom;
256+
257+
// COPYFROMTO(Var.LoopMom[LoopIndex], CurrMom);
258+
CurrMom = Var.LoopMom[LoopIndex];
259+
263260
if (LoopIndex == 0) {
264261
Prop = ShiftExtK(Var.CurrExtMomBin, NewExtMomBin);
265-
COPYFROMTO(Var.ExtMomTable[NewExtMomBin], Var.LoopMom[LoopIndex]);
262+
// COPYFROMTO(Var.ExtMomTable[NewExtMomBin], Var.LoopMom[LoopIndex]);
263+
Var.LoopMom[LoopIndex] = Var.ExtMomTable[NewExtMomBin];
266264
} else {
267265
Prop = ShiftK(CurrMom, Var.LoopMom[LoopIndex]);
268266
}
269-
if (LoopIndex == 0 && norm2(Var.LoopMom[LoopIndex]) > Para.MaxExtMom) {
270-
COPYFROMTO(CurrMom, Var.LoopMom[LoopIndex]);
267+
if (LoopIndex == 0 && Var.LoopMom[LoopIndex].norm() > Para.MaxExtMom) {
268+
Var.LoopMom[LoopIndex] = CurrMom;
269+
// COPYFROMTO(CurrMom, Var.LoopMom[LoopIndex]);
271270
return;
272271
}
273272

@@ -280,7 +279,8 @@ void markov::ChangeMomentum() {
280279
if (LoopIndex == 0)
281280
Var.CurrExtMomBin = NewExtMomBin;
282281
} else {
283-
COPYFROMTO(CurrMom, Var.LoopMom[LoopIndex]);
282+
Var.LoopMom[LoopIndex] = CurrMom;
283+
// COPYFROMTO(CurrMom, Var.LoopMom[LoopIndex]);
284284
Weight.RejectChange(*Var.CurrGroup);
285285
}
286286
};
@@ -332,7 +332,7 @@ double markov::RemoveOldK(momentum &OldMom) {
332332
double dK = Para.Kf / sqrt(Para.Beta) / 4.0;
333333
if (dK > Para.Kf / 2)
334334
dK = Para.Kf / 2; // to avoid dK>Kf situation
335-
double KAmp = norm2(OldMom);
335+
double KAmp = OldMom.norm();
336336
if (KAmp < Para.Kf - dK || KAmp > Para.Kf + dK)
337337
// Kf-dK<KAmp<Kf+dK
338338
return 0.0;
@@ -357,13 +357,14 @@ double markov::ShiftK(const momentum &OldMom, momentum &NewMom) {
357357
double x = Random.urn();
358358
double Prop;
359359
if (x < 1.0 / 3) {
360-
COPYFROMTO(OldMom, NewMom);
360+
// COPYFROMTO(OldMom, NewMom);
361+
NewMom = OldMom;
361362
int dir = Random.irn(0, D - 1);
362363
double STEP = Para.Beta > 1.0 ? Para.Kf / Para.Beta * 3.0 : Para.Kf;
363364
NewMom[dir] += STEP * (Random.urn() - 0.5);
364365
Prop = 1.0;
365366
} else if (x < 2.0 / 3) {
366-
double k = norm2(OldMom);
367+
double k = OldMom.norm();
367368
if (k < EPS)
368369
Prop = 0.0;
369370

‎src/markov.h

+4-3
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include "utility/rng.h"
66
#include "weight.h"
77
#include <string>
8+
#include <unordered_map>
89
#include <vector>
910

1011
namespace mc {
@@ -36,11 +37,11 @@ class markov {
3637
int DynamicTest();
3738

3839
private:
39-
// Observables
4040
// polarizatoin for each group
41-
vector<polar> Polar;
41+
unordered_map<int, polar> Polar;
42+
4243
// polarizatoin for each group at the zero momentumr;
43-
vector<double> PolarStatic;
44+
unordered_map<int, double> PolarStatic;
4445

4546
// MC updates
4647

0 commit comments

Comments
 (0)