Class: CodeRunner::Corfu::Optimisation

Inherits:
Object
  • Object
show all
Includes:
GSL::MultiMin
Defined in:
lib/corfucrmod/optimisation.rb

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(optimised_quantity, optimisation_spec) ⇒ Optimisation

Returns a new instance of Optimisation.



29
30
31
32
33
34
35
36
37
38
39
40
41
 
# File 'lib/corfucrmod/optimisation.rb', line 29

def initialize(optimised_quantity, optimisation_spec)
  #@folder = folder
  @optimised_quantity = optimised_quantity
  @optimisation_spec = optimisation_spec
  @optimisation_variables = optimisation_spec.map{|code, hash| hash.map{|var, pars| [code, var]}}.flatten(1)
  @optimisation_starts    = optimisation_spec.map{|code, hash| hash.map{|var, pars| pars[0]}}.flatten(1)
  @optimisation_steps     = optimisation_spec.map{|code, hash| hash.map{|var, pars| pars[1]}}.flatten(1)
  #@runner = CodeRunner.fetch_runner(
  #p ['optimisation_variables', @optimisation_variables]
  @results_hash = {}
  @first_call = true
  @ifspppl_converged = false
end

Instance Attribute Details

#gs_runner ⇒ Object

Returns the value of attribute gs_runner.



26
27
28
 
# File 'lib/corfucrmod/optimisation.rb', line 26

def gs_runner
  @gs_runner
end

#optimisation_spec ⇒ Object (readonly)

optimisation_spec is a hash of => {:variable => [initial_guess, dimension_scale_factor]} dimension_scale_factor is just some estimate of the length scale in which the result varies significantly code_name is either trinity or chease (both can be used simultaneously)



23
24
25
 
# File 'lib/corfucrmod/optimisation.rb', line 23

def optimisation_spec
  @optimisation_spec
end

#optimisation_variables ⇒ Object (readonly)

Returns the value of attribute optimisation_variables.



24
25
26
 
# File 'lib/corfucrmod/optimisation.rb', line 24

def optimisation_variables
  @optimisation_variables
end

#parameters_obj ⇒ Object

Returns the value of attribute parameters_obj.



27
28
29
 
# File 'lib/corfucrmod/optimisation.rb', line 27

def parameters_obj
  @parameters_obj
end

#results_hash ⇒ Object

Returns the value of attribute results_hash.



28
29
30
 
# File 'lib/corfucrmod/optimisation.rb', line 28

def results_hash
  @results_hash
end

#trinity_runner ⇒ Object

Returns the value of attribute trinity_runner.



25
26
27
 
# File 'lib/corfucrmod/optimisation.rb', line 25

def trinity_runner
  @trinity_runner
end

Instance Method Details

#dimension ⇒ Object 



42
43
44
 
# File 'lib/corfucrmod/optimisation.rb', line 42

def dimension
  @optimisation_variables.size
end

#func(v) ⇒ Object 



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
 
# File 'lib/corfucrmod/optimisation.rb', line 79

def func(v)
  val = nil
  count = 1
  val_old, repeat = func_actual(v)
  loop do
    val, repeat, trinity_is_converged = func_actual(v)
    if trinity_is_converged # This means that Trinity has reached steady state

      # Now we check if the loop over recalculating the GS equation
      # has converged
      if ((val_old - val)/val).abs < @parameters_obj.convergence
        if @parameters_obj.use_ifspppl_first and not @ifspppl_converged
          @ifspppl_converged = true
        else
          break
        end
      end
      break if count > @parameters_obj.max_func_evals

      # If geometry is being evolved internally in Trinity, we don't
      # need to run this loop
      break if not repeat
      val_old = val
      #break if not repeat or (not @first_call and count > 1) or count > 4
      count += 1 unless @parameters_obj.use_ifspppl_first and not @ifspppl_converged
    end
  end
  @first_call = false
  return val
end

#func_actual(v) ⇒ Object 



109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
 
# File 'lib/corfucrmod/optimisation.rb', line 109

def func_actual(v)
  eputs 'Starting func'
  @id||=1 # Id of the current trinity run which should always be equal to the number of func_actual calls
  repeat = false
  print_pars = {}
  pars = {}
  pars[:gs] = {}
  pars[:trinity] = {}
  pars[:trinity].absorb(@parameters_obj.trinity_pars) if @parameters_obj.trinity_pars
  pars[:gs].absorb(@parameters_obj.gs_pars) if @parameters_obj.gs_pars
  for i in 0...v.size
    code, varname = @optimisation_variables[i]
    val = v[i]
    code = case code
           when :chease, :ecom
             :gs
           else
             :trinity
           end
    pars[code][varname] = val
    print_pars[code] ||={}
    print_pars[code][varname] = val
  end
  if not @first_run_done
    pars[:trinity][:ntstep] = @parameters_obj.ntstep_first
    #pars[:trinity][:nifspppl_initial] = -1
    #pars[:trinity][:niter] = 2

    # The line below assumes that the whole first run is done
    # with ifspppl, in which case we don't want the timestep to get
    # too large for the gryfx run afterwards.
    #pars[:trinity][:ntdelt_max] = 0.05
    #pars[:trinity][:convergetol] = -1.0
  else
    pars[:trinity][:ntstep] = @parameters_obj.ntstep
    #pars[:trinity][:nifspppl_initial] = -1
    #pars[:trinity][:niter] = 3
    #pars[:trinity][:ntdelt_max] = 10.0
  end
  if @parameters_obj.use_ifspppl_first and not @ifspppl_converged
    pars[:trinity][:nifspppl_initial] = pars[:trinity][:ntstep]
  else
    pars[:trinity][:nifspppl_initial] = -1
  end


  # Must fix this soon!
  if @parameters_obj.gs_code == 'chease'
    #pars[:gs][:ap] = [0.3,0.5,0.4,0.0,0.4,0.0,0.0]
    #pars[:gs][:at] = [0.16,1.0,1.0,-1.1,-1.1]
  end

  trinity_runner.run_class.instance_variable_set(:@mpi_communicator, MPI::Comm::WORLD)
  #if false and trinity_runner.run_list.size > 0
  #else
  if @first_run_done
    #crun.nppfun=4
    #crun.neqdsk=0
    #crun.expeq_file = trinity_runner.run_list[@id]
  end
  if not @replay
    if not @first_run_done and gs_runner.run_list.size > 0
      #This means we are in a restart
      @replay = true
      @nrun = 1
      if @parameters_obj.delete_final_run
        eputs 'Removing final run: ' + trinity_runner.run_list.keys.max.to_s
        trinity_runner.conditions =  'id==' + trinity_runner.run_list.keys.max.to_s
        trinity_runner.destroy no_confirm: true
        gs_runner.conditions =  'id==' + gs_runner.run_list.keys.max.to_s
        gs_runner.destroy no_confirm: true
      end
    else
      @nrun = nil
      @replay = false
    end
  end
  if @replay
    eputs 'Replaying: ' + @nrun.to_s
    run = trinity_runner.run_list[@nrun]
    if not run
      eputs 'Ending replay at ' + @nrun.to_s
      @replay = false
      @id = @gsid = @nrun-1
    else
      @id = @gsid = @nrun
    end
    @nrun += 1
  end
  Hash.phoenix('func_calls.rb') do |hash|
    hash[@id+1] ||={}
    hash[@id+1][:variables] = v.dup
    hash[@id+1][:print_pars] = print_pars
    hash
  end
  #raise "WORK!!"
  eputs "Written parameters"
  if not @replay
    eputs "Not replaying... starting GS and Trinity runs"
    remaining_wall_mins = (
      @parameters_obj.wall_mins - @parameters_obj.wall_mins_margin -
      (Time.now.to_i - @parameters_obj.start_time).to_f / 60.0
    )
    eputs "Remaining wall mins #{remaining_wall_mins}, wall mins #{@parameters_obj.wall_mins}, start time #{@parameters_obj.start_time}, time #{Time.now.to_i}, margin #{@parameters_obj.wall_mins_margin}"
    if remaining_wall_mins < @parameters_obj.wall_mins_margin
      eputs "Run out of time"
      throw(:out_of_time)
    end
    eputs "Starting real run, @id = ",@id
    # Create and initialize the gs run
    gsrun = gs_runner.run_class.new(gs_runner)
    raise "No gs_defaults strings" unless @parameters_obj.gs_defaults_strings.size > 0
    @parameters_obj.gs_defaults_strings.each{|prc| gsrun.instance_eval(prc)}
    if gs_runner.run_list.size > 0
      #gsrun.restart_id = @gsid
      if @parameters_obj.gs_code == 'chease'
        last_converged = @id
        #We need to find the last converged Trinity run to use as the pressure profile.
        until last_converged == 0 or trinity_runner.combined_run_list[last_converged].is_converged?
          eputs "Run #{last_converged} not converged"
          last_converged -= 1
        end
        #unless (@parameters_obj.use_previous_pressure==0 and not @first_trinity_run_completed)
        unless last_converged == 0
          eputs "Using previous pressure profile"
          # We give CHEASE the pressure profile from the previous run.
          pars[:gs][:nppfun] = 4
          pars[:gs][:nfunc] = 4
          #if prid = @parameters_obj.use_previous_pressure and not @first_trinity_run_completed
          # If the last trinity run did not converge we may want to run exactly
          # the same case again, and in particular use the pressure profile from
          # the previous Trinity  run as input (as an unconverged pressure profile
          # can lead to a wacky GS solution)
          #gsrun.expeq_in=trinity_runner.combined_run_list[prid].directory + '/chease/EXPEQ.NOSURF'
          #else
          gsrun.expeq_in=trinity_runner.combined_run_list[last_converged].directory + '/chease/EXPEQ.NOSURF'
          #end
          # Don't optimise presssure profile.
          pars[:gs][:nblopt] = 0
        end
      end
    end
    gsrun.update_submission_parameters(pars[:gs].inspect, false)

    # Create and initialize the trinity run
    run = trinity_runner.run_class.new(trinity_runner)
    raise "No trinity_defaults_strings" unless @parameters_obj.trinity_defaults_strings.size > 0
    run.instance_variable_set(:@set_flux_defaults_procs, []) unless run.instance_variable_get(:@set_flux_defaults_procs)
    @parameters_obj.trinity_defaults_strings.each{|prc| run.instance_eval(prc)}
    run.update_submission_parameters(pars[:trinity].inspect, false)

    #if @parameters_obj.gs_code == 'chease' and (run.evolve_geometry and run.evolve_geometry.fortran_true?)
    #pars[:gs][:nblopt] = 0
    #end
    gs_runner.submit(gsrun)
    gsrun = gs_runner.run_list[@gsid = gs_runner.max_id]
    gsrun.recheck
    gs_runner.update
    #gs_runner.print_out(0)
    #FileUtils.cp(gsrun.directory + '/ogyropsi.dat', trinity_runner.root_folder + '/.')

    run.gs_folder = gsrun.directory
    while (
      not FileTest.exist? run.gs_folder + '/ogyropsi.dat' or
      File.read(run.gs_folder + '/ogyropsi.dat') =~ /nan/i
    )
      #eputs "GS solver failed: using previous solution"

      gs_runner.conditions = 'id == ' + @gsid.to_s
      gs_runner.destroy(no_confirm: true)
      gs_runner.conditions = nil
      eputs "GS solver failed for #{v.inspect}: returning 10000"
      return [10000, false]
      #run.gs_folder = gs_runner.run_list[@gsid -= 1].directory
    end
    #run.evolve_geometry = ".true."
    eputs ['Set gs_folder', run.gs_folder]
    trinity_runner.run_class.instance_variable_set(:@delay_execution, true)
    if trinity_runner.run_list.size > 0
      if (old_run = trinity_runner.run_list[@id]).is_converged?
        eputs "Previous trinity run #@id converged: using profile as initial condition"
        run.init_option="trinity"
        ## We have to copy the file because we need to be able to access the original from R
        #FileUtils.cp("../id_#@id/#{old_run.run_name}.out.nc", "../id_#@id/#{old_run.run_name}_copy.out.nc")
        run.init_file = "../id_#@id/#{old_run.run_name}.out.nc"
        run.init_time = old_run.new_netcdf_file.var("t").get[-1]
        old_run.new_ncclose rescue nil # Make sure the file is closed
      else
        eputs "Previous trinity run #@id not converged: restarting"
        run.restart_id = @id
      end 
    end
    eputs 'Submitting run'
    run.wall_mins = remaining_wall_mins
    trinity_runner.submit(run)
    run = trinity_runner.run_list[@id = trinity_runner.max_id]


    comm = MPI::Comm::WORLD
    arr = NArray.int(1)
    arr[0] = 1
    eputs 'Sending message'
    comm.Bcast(arr,0)
    comm.Bcast_string(run.directory)
    comm.Bcast_string(run.run_name)
    eputs 'Running trinity'
    Dir.chdir(run.directory){run.run_trinity(run.run_name+'.trin', comm)}
    eputs 'Rechecking'
    trinity_runner.update
    eputs 'Queue', run.queue_status

    trinity_runner.update
    @first_trinity_run_completed
  end # if not @replay
  #trinity_runner.print_out(0)
  Dir.chdir(run.directory) do
    run.recheck
    run.status = :Complete
    run.get_global_results
    if (run.evolve_geometry and run.evolve_geometry.fortran_true?)
      repeat = false
    else
      repeat = true
    end
  end
  result =  run.instance_eval(@optimised_quantity)
  p ['result is ', result, 'repeat: ', repeat]
  @first_run_done = true
  @results_hash[:func_calls] ||=[]
  @results_hash[:func_calls].push [print_pars, result]
  Hash.phoenix('func_calls.rb') do |hash|
    hash[@id] ||={}
    hash[@id][:result] = result
    hash[@id][:repeat] = repeat
    hash[@id][:is_converged] = run.is_converged?
    hash
  end
  return [-result, repeat, run.is_converged?]
  #end


  #v.square.sum
end

#serial_optimise(optimisation_method, parameters_obj) ⇒ Object 



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
 
# File 'lib/corfucrmod/optimisation.rb', line 45

def serial_optimise(optimisation_method, parameters_obj)
  optimisation_meth = case optimisation_method
                      when :simplex
                        FMinimizer::NMSIMPLEX
                      else
                        raise "Unknown optimisation_method"
                      end
  @results_hash[:start_time] = Time.now.to_i
  opt = FMinimizer.alloc(optimisation_meth, @optimisation_variables.size)
  @parameters_obj = parameters_obj
  catch(:out_of_time) do
    func = Proc.new{|v, optimiser| optimiser.func(v)}
    eputs 'Created func'
    gsl_func = Function.alloc(func, dimension)
    eputs 'Allocated gsl_func'
    gsl_func.set_params(self)
    eputs 'Set params'
    opt.set(gsl_func, @optimisation_starts.to_gslv, @optimisation_steps.to_gslv)
    eputs 'Set func and starting iteration'
    parameters_obj.nit.times do |i|
      opt.iterate
      p ['status', opt.x, opt.minimum, i, parameters_obj.nit]
      @results_hash[:iterations] ||= []
      @results_hash[:iterations].push [opt.x.dup, opt.minimum]
      @results_hash[:elapse_mins] = (Time.now.to_i - @results_hash[:start_time]).to_f/60
      @results_hash[:current_time] = Time.now.to_i
      File.open('results', 'w'){|f| f.puts @results_hash.pretty_inspect}
    end
    eputs 'Optimisation complete'
  end
  eputs 'Optimisation ended'
  #MPI.Finalize

end