forked from DSSAT/glue
-
Notifications
You must be signed in to change notification settings - Fork 0
/
OptimalParameterSet.r
122 lines (99 loc) · 4.64 KB
/
OptimalParameterSet.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
##This is the function to get the optimal parameter set for future research.
OptimalParameterSet<-function (GLUEFlag, OD, DSSATD, CropName, CultivarID, CultivarName, GenotypeFileName, TotalParameterNumber)
{
##1.Select the optimal parameter set.
if(GLUEFlag==1 | GLUEFlag==3)
{
eval(parse(text=paste('Optimal<-read.table("',OD,
'/PosteriorDistribution_2.txt", header=TRUE,comment.char="")',sep = '')));
OptimalParameterSet<-Optimal["MaxProbability",];
} else
{
eval(parse(text=paste('Optimal<-read.table("',OD,
'/PosteriorDistribution_1.txt", header=TRUE,comment.char="")',sep = '')));
OptimalParameterSet<-Optimal["MaxProbability",];
}
##2.Change partial cultivar file above with the values of optimal parameter set.
## This cultivar then can be chosen by model users to do future simulations. It requires model
## users to copy this cultivar to replace the corresponding cultivar in the genotype file under "DSSAT048/Genotype".
eval(parse(text=paste('GenotypeFilePath="',OD,'/',GenotypeFileName,'.CUL"',sep = '')));
GenotypeFile<-readLines(GenotypeFilePath, n=-1);
# Read the genotype file in the working directory.
CultivarID<-paste("^",CultivarID, sep='');
CultivarAddress<-grep(pattern=CultivarID, GenotypeFile);
# Read the line that has the information of the cultivar that is under estimation.
OldLine<-GenotypeFile[CultivarAddress];
#Get the line according to the line number.
if(CropName != "SC")
{
Step<-6;
ValuePosition1<-(38-Step);
ValuePosition2<-(42-Step);
#Set the starting and end points of parameter locations.
for (i in 1:TotalParameterNumber)
{
ValuePosition1<-ValuePosition1+Step;
ValuePosition2<-ValuePosition2+Step;
eval(parse(text = paste("Parameter<-OptimalParameterSet[1,",i,"]",sep = '')));
#To solve the format problem for parameters with negative values. Modified by He, 2015-6-18.
if(Parameter < 0 & Parameter > -1.0) #
{ #
ParameterFormat<-sprintf('%1.3f', Parameter); #
ParameterFormat<-paste(substring(ParameterFormat,1,1), substring(ParameterFormat,3), sep=''); #
} else if (Parameter <= -1.0 & Parameter > -10.0) #
{ #
ParameterFormat<-sprintf('%2.2f', Parameter); #
} else if (Parameter <= -10.0 & Parameter > -100.0) #
{ #
ParameterFormat<-sprintf('%3.1f', Parameter); #
} #
if(Parameter>=0 & Parameter<10)
{
ParameterFormat<-sprintf('%1.3f', Parameter);
} else if (Parameter>=10 & Parameter<100)
{
ParameterFormat<-sprintf('%2.2f', Parameter);
} else if (Parameter >= 100)
{
ParameterFormat<-sprintf('%3.1f', Parameter);
}
substr(OldLine, ValuePosition1, ValuePosition2)<-ParameterFormat;
}
} else
{
Step<-15;
#chp modified
ValuePosition1<-(47-Step);
ValuePosition2<-(61-Step);
#Set the starting and end points of parameter locations.
for (i in 1:TotalParameterNumber)
{
ValuePosition1<-ValuePosition1+Step;
ValuePosition2<-ValuePosition2+Step;
eval(parse(text = paste("Parameter<-OptimalParameterSet[1,",i,"]",sep = '')));
if(Parameter>=0 & Parameter<10)
{
ParameterFormat<-sprintf('%1.3f', Parameter);
} else if (Parameter>=10 & Parameter<100)
{
ParameterFormat<-sprintf('%2.2f', Parameter);
# chp added extra format statement for values between 100 and 1000
} else if (Parameter>=100 & Parameter<1000)
{
ParameterFormat<-sprintf('%3.1f', Parameter);
} else
{
ParameterFormat<-sprintf('%4.0f', Parameter);
}
#print(ParameterFormat);
substr(OldLine, ValuePosition1, ValuePosition2)<-' ';# Delete initial values.
substr(OldLine, ValuePosition1, ValuePosition2)<-ParameterFormat;
}
}
NumberOfIllegalCharacters<-grep(pattern="/", CultivarName);
if(length(NumberOfIllegalCharacters)!=0) CultivarName<-gsub(pattern="/", "_", CultivarName);
##For some cultivars, they have illegal characters such as /, *, # in theri name. This should be changed.
eval(parse(text=paste("NewGenotypeFilePath='",OD,"/",CropName,"",CultivarName,".CUL'",sep = '')));
write(OldLine, file=NewGenotypeFilePath);
#Save the new genotype file as "cul" file in the DSSAT directory.
}