plot.jl (6853B)
1 #!/usr/bin/env julia 2 import ArgParse 3 ENV["MPLBACKEND"] = "Agg" 4 import PyPlot 5 import LsqFit 6 7 verbose = false 8 9 function parse_command_line() 10 s = ArgParse.ArgParseSettings() 11 ArgParse.@add_arg_table s begin 12 "--nruns" 13 help = "number of runs in ensemble" 14 arg_type = Int 15 default = 1 16 "id" 17 help = "simulation id" 18 required = true 19 end 20 return ArgParse.parse_args(s) 21 end 22 23 function report_args(parsed_args) 24 println("Parsed args:") 25 for (arg,val) in parsed_args 26 println(" $arg => $val") 27 end 28 end 29 30 function main() 31 parsed_args = parse_command_line() 32 #report_args(parsed_args) 33 34 sim_id = parsed_args["id"] 35 #info(sim_id) 36 nruns = parsed_args["nruns"] 37 t = Vector 38 #t_jam = ones(nruns)*parsed_args["total_hours"]*60.*60.*2. 39 t_jam = ones(nruns)*1e16 40 nskipped = 0 41 iteration_skipped = zeros(Bool, nruns) 42 avg_coordination_number = Vector[] 43 44 PyPlot.figure(figsize=(5,4)) 45 46 for i=1:nruns 47 seed = i 48 filename = parsed_args["id"] * "-seed" * string(seed) * "-ice-flux.txt" 49 #info(filename) 50 if !isfile(filename) 51 nskipped += 1 52 iteration_skipped[i] = true 53 push!(avg_coordination_number, []) 54 continue 55 end 56 indata = readdlm(filename) 57 #t, ice_flux = indata[1,:], indata[2,:] 58 t, ice_flux, avg_coordination_no = 59 indata[1,10:end], indata[2,10:end], indata[3,10:end] # skip first 10 vals 60 t -= t[1] 61 push!(avg_coordination_number, avg_coordination_no) 62 PyPlot.plot(t/(60.*60.), ice_flux) 63 64 time_elapsed_while_jammed = 0. 65 for it=(length(t) - 1):-1:1 66 if ice_flux[it] ≈ ice_flux[end] 67 time_elapsed_while_jammed = t[end] - t[it] 68 end 69 end 70 71 if time_elapsed_while_jammed > 60.*60. 72 t_jam[i] = t[end] - time_elapsed_while_jammed 73 #info("simulation $i jammed at t = $(t_jam[i]/(60.*60.)) h") 74 end 75 76 end 77 #PyPlot.title(parsed_args["id"]) 78 PyPlot.xlabel("Time [h]") 79 PyPlot.ylabel("Cumulative ice throughput [kg]") 80 PyPlot.tight_layout() 81 PyPlot.title("(a)") 82 83 PyPlot.savefig(parsed_args["id"] * ".png") 84 PyPlot.savefig(parsed_args["id"]) 85 86 PyPlot.clf() 87 jam_fraction = zeros(length(t)) 88 ensemble_avg_coordination_number = zeros(length(t)) 89 for it=1:length(t) 90 for i=1:nruns 91 if iteration_skipped[i] 92 continue 93 end 94 if t_jam[i] <= t[it] 95 jam_fraction[it] += 1./float(nruns - nskipped) 96 end 97 if it > length(avg_coordination_number[i]) 98 continue 99 end 100 ensemble_avg_coordination_number[it] += 101 avg_coordination_number[i][it]/float(nruns - nskipped) 102 end 103 end 104 survived_fraction = 1. - jam_fraction 105 t_hours = t/(60.*60.) 106 107 # fit function of exponential decay to survived_fraction 108 survival_model(t_hours, tau) = exp.(-t_hours/tau[1]) 109 tau0 = [2.] # initial guess of tau [hours] 110 I = find(jam_fraction) # find first where survived_fraction > 0. 111 first_idx = 1 112 if length(I) > 0 113 first_idx = I[1] 114 end 115 if t_hours[first_idx] > 6. 116 first_idx = 1 117 end 118 t_hours_trim = linspace(0., t_hours[end], 119 length(survived_fraction[first_idx:end])) 120 fit = LsqFit.curve_fit(survival_model, 121 t_hours_trim, 122 survived_fraction[first_idx:end], 123 tau0) 124 125 covar = LsqFit.estimate_covar(fit) 126 std_error = sqrt.(abs.(diag(covar))) 127 sample_std_dev = std_error .* sqrt(length(t_hours_trim)) 128 errors = LsqFit.estimate_errors(fit, 0.95) 129 130 if fit.converged 131 #print_with_color(:green, 132 #"tau = $(fit.param[1]) h, 95% s = $(sample_std_dev[1]) h\n") 133 #"tau = $(fit.param[1]) h, 95% CI = $(errors[1]) h\n") 134 println("$(sim_id) $(fit.param[1]) $(sample_std_dev[1])") 135 #label = @sprintf("\$P_s(t, \\tau = %4.3g\$ h), \$SE_\\tau = %4.3g\$ h", 136 #fit.param[1], std_error[1]) 137 label = @sprintf("\$P_s(t, T = %4.3g\$ h), \$s_T = %4.3g\$ h", 138 fit.param[1], sample_std_dev[1]) 139 140 # 95% CI range for normal distribution 141 #lower_CI = survival_model(t_hours, fit.param - 1.96*std_error[1]) 142 #upper_CI = survival_model(t_hours, fit.param + 1.96*std_error[1]) 143 144 # +/- one sample standard deviation 145 lower_SD = survival_model(t_hours, fit.param - sample_std_dev[1]) 146 upper_SD = survival_model(t_hours, fit.param + sample_std_dev[1]) 147 148 #PyPlot.plot(t_hours, survival_model(t_hours, fit.param - errors[1]), 149 #":C1", linewidth=1) 150 #PyPlot.plot(t_hours, survival_model(t_hours, fit.param + errors[1]), 151 #":C1", linewidth=1) 152 153 PyPlot.fill_between(t_hours + t_hours[first_idx], 154 lower_SD, upper_SD, facecolor="C1", 155 interpolate=true, alpha=0.33) 156 157 PyPlot.plot(t_hours + t_hours[first_idx], 158 survival_model(t_hours, fit.param), 159 "--C1", linewidth=2, label=label) 160 161 else 162 warn(parsed_args["id"] * ": LsqFit.curve_fit did not converge") 163 end 164 165 # Plot DEM 166 PyPlot.plot(t_hours, survived_fraction, 167 "-C0", linewidth=2, 168 label="DEM (\$n = $(nruns - nskipped) \$)") 169 170 if fit.param[1] > 30. 171 PyPlot.legend(loc="lower right", fancybox="true") 172 else 173 PyPlot.legend(loc="upper right", fancybox="true") 174 end 175 176 #PyPlot.title(parsed_args["id"] * ", N = " * string(nruns - nskipped)) 177 PyPlot.title("(b)") 178 PyPlot.xlabel("Time [h]") 179 PyPlot.ylabel("Probability of survival \$P_s\$ [-]") 180 PyPlot.xlim([minimum(t_hours), maximum(t_hours)]) 181 PyPlot.ylim([-0.02, 1.02]) 182 #PyPlot.ylim([0., 1.]) 183 184 PyPlot.tight_layout() 185 PyPlot.savefig(parsed_args["id"] * "-survived_fraction.pdf") 186 PyPlot.savefig(parsed_args["id"] * "-survived_fraction.pdf.png") 187 writedlm(parsed_args["id"] * "-survived_fraction.txt", 188 [t_hours, survived_fraction, ensemble_avg_coordination_number]) 189 190 # Plot ensemble-mean coordination number 191 PyPlot.clf() 192 PyPlot.plot(t_hours, ensemble_avg_coordination_number, 193 "-C0", linewidth=2) 194 PyPlot.xlabel("Time [h]") 195 PyPlot.ylabel("Avg. coordination number \$Z\$ [-]") 196 PyPlot.xlim([minimum(t_hours), maximum(t_hours)]) 197 PyPlot.ylim([-0.02, 5.02]) 198 PyPlot.tight_layout() 199 PyPlot.savefig(parsed_args["id"] * "-coordination_number.pdf") 200 PyPlot.savefig(parsed_args["id"] * "-coordination_number.pdf.png") 201 202 end 203 204 main()