[Ifeffit] Hard restraints?

[Matt Newville]

Hi Scott, Theanne,

> .... We can't just put an abs() around the def'd values,
> because that means our constraints based on stoichiometry and
> the relative number of sites are being applied incorrectly. A
> restraint doesn't really seem appropriate either.
>
> Does anyone have an idea of how to deal with this?

In a simplified case, you'd like to constrain two coordination
numbers to be positive _and_ to add up to some value, right?  
  guess n1_var = 5
  def   n1     = abs(n1_var)
  def   n2     = 12 - n1

or, if you're worried about n1 > 12, n2<0, you could do: 
  guess n1_var = 5
  def   n1     = max(0,min(n1_var,12))
  def   n2     = 12 - n1

Is that good enough, or am I not seeing the whole complicated
picture?

> If no one does, perhaps there should be an ifeffit option for
> a "hard restraint." This would be something like the max and
> min functions, except that it would operate by putting a huge
> penalty to the chi-square when outside the range. Just to be
> clear, it would work something like this: y has been def'd to
> x + z. Then y is also given a hard restraint that it cannot be
> less than 0. The fit would then generate a positive value of y
> that was not less than 0 by varying x and z in such a way that
> it is still true that y = x + z.

I think this is possible without a change to the code.  You can
add a restraint penalty when some value is outside a range
[lo_val,hi_val].  It's a bit scary looking, but this will do it:

  lo_val   = 0.7    ## lower bound
  up_val   = 1.0    ## upper bound
  scale    = 1.000  ## scale can be increased to weight this restraint
  target   = 0.8    ## this could be a variable, say S02

  ## use penalty as a restraint in feffit()
  def  penalty  = (max((target-lo_val), abs(target-lo_val)) - 
                   min((target-up_val),-abs(target-up_val)) + 
                   lo_val - up_val)*scale

This penalty can be used as a restraint, and will be 0 when the 
target value is between the lower and upper bound.

To set a penalty only for going below the lower bound, it's just
  def penalty = max((target-lo_val),abs(target-lo_val))-(target- lo_val)

and if the lower bound is 0, it's even simpler:
  def penalty = max(target,abs(target))-target

Even that is scary enough, and the whole thing may suggest a
change to the code to make a simple function that provides a
restraint penalty when a variable is outside some user-selected
bounds, perhaps as
    bound(x, lo_val, hi_val)  # not implemented!

Is this way of setting a penalty what you had in mind, or am I
missing the point?


--Matt

PS: 

##########
# To visualize this restraint, an array of penalties can be 
# plotted v. an array of targe values:
  lo_val   = 0.7
  up_val   = 1.0
  scale    = 1.000

  # array of target values
  m.vals = range(-1,3,0.05)

  m.lo_bound = scale*(m.vals-up_val)
  m.up_bound = scale*(m.vals-lo_val)

  m.lo_penalty =  max(m.up_bound,  abs(m.up_bound)) - m.up_bound
  m.up_penalty = -min(m.lo_bound, -abs(m.lo_bound)) + m.lo_bound

  m.penalty =  (max(m.up_bound,  abs(m.up_bound)) 
              - min(m.lo_bound, -abs(m.lo_bound)) + scale*(lo_val - up_val))

  newplot m.vals, m.penalty, style=linespoints2
  plot m.vals, m.lo_penalty

##########