Inversion example

The script Interp.py does out-of-core interpolation. The program Interp.x interpolates from an irregular to a regular mesh. The program reg.x convolves with a three point filter. We then can turn these program into operators.

interp_op=SEP.op_oc_serial.operator("Interp.x")
reg      =SEP.op_oc_serial.operator("reg.x")

We then extend the solver object. We require the operator op and regularization operator reg. (n1,o1,d1).

class solver(SEP.solv_reg.solver):
  def __init__(self,op,reg):
    SEP.solv_reg.solver.__init__("SOLVER")
    self.add_param("op",op)
    self.add_param("reg",reg)

We require the user to specify the model model, the data data, $\epsilon$ eps, and the number of iterations niter, and the dimensions of the output space

  def add_doc_params(self):
    self.add_doc_param("eps",1.,"Epsilon")
    self.add_doc_param("niter",1,"Number of iterations")
    self.add_doc_param("model",doc="Model (output)")
    self.add_doc_param("data",doc="Data (input)")
    self.add_doc_param("n1",doc="number of samples (axis 1)")
    self.add_doc_param("o1",doc="First sample (axis 1)")
    self.add_doc_param("d1",doc="Sampling  (axis 1)")

We create the data vector from a file and create the model vector.

  def prep_run(self,restart=None):
    self.add_param("data",SEP.vec_sfloat.vector(tag="data"))
    self.add_param("model",SEP.vec_sfloat.vector(name="model"))
    self.param("model").set_axis(1, self.param("n1").
      self.param("o1",self.param("d1")))
    self.param("model").zero()

We add the domain and range vectors for the operator and the regularization operator.

    self.param("op").add_param("domain",self.main.param("model"))
    self.param("op").add_param("range",self.main.param("data"))
    self.param("reg").add_param("domain",self.main.param("model"))
    self.param("reg").add_param("range",self.main.param("model"))

Finally we create the step operator and run the prep_run function for the operators and the regularization solver.

    self.add_param("cgstep",SEP.cgstep.cgstep("cgstep"))
    self.param("op").prep_run(restart)
    self.param("reg").prep_run(restart)
    SEP.solv_reg.solver.prep_run(self,restart)

We create the solver object, the program object, read the parameters, and run the job.

solv=solver(interp_op,reg_op)
                                                                                    
program=SEP.opt_prog.prog("Interpolate  ",
   "Interp.py model= data= ",
    [interp_op,reg_op,solv],
     ["Interpolatation  using the SEP python library"])
program.get_options()
program.prep_run()