def impose_reweighted_variance(v, samples, weights=None, solver=None):
    """impose a variance on a list of points by reweighting weights"""
    ndim = len(samples)
    if weights is None:
        weights = [1.0/ndim] * ndim
    if solver is None or solver == 'fmin':
        from mystic.solvers import fmin as solver
    elif solver == 'fmin_powell':
        from mystic.solvers import fmin_powell as solver
    elif solver == 'diffev':
        from mystic.solvers import diffev as solver
    elif solver == 'diffev2':
        from mystic.solvers import diffev2 as solver
    norm = sum(weights)
    m = mean(samples, weights)

    inequality = ""
    equality = ""; equality2 = ""; equality3 = ""
    for i in range(ndim):
        inequality += "x%s >= 0.0\n" % (i) # positive
        equality += "x%s + " % (i)         # normalized
        equality2 += "%s * x%s + " % (float(samples[i]),(i)) # mean
        equality3 += "x%s*(%s-%s)**2 + " % ((i),float(samples[i]),m) # var

    equality += "0.0 = %s\n" % float(norm)
    equality += equality2 + "0.0 = %s*%s\n" % (float(norm),m)
    equality += equality3 + "0.0 = %s*%s\n" % (float(norm),v)

    penalties = generate_penalty(generate_conditions(inequality))
    constrain = generate_constraint(generate_solvers(solve(equality)))

    def cost(x): return sum(x)

    results = solver(cost, weights, constraints=constrain, \
                     penalty=penalties, disp=False, full_output=True)
    wts = list(results[0])
    _norm = results[1] # should have _norm == norm
    warn = results[4]  # nonzero if didn't converge

    #XXX: better to fail immediately if xlo < m < xhi... or the below?
    if warn or not almostEqual(_norm, norm):
        print "Warning: could not impose mean through reweighting"
        return None #impose_variance(v, samples, weights), weights

    return wts #samples, wts  # "mean-preserving"

    return -(sin(2*pi*x0)**3 * sin(2*pi*x1)) / (x0**3 * (x0 + x1))

bounds = [(0,10)]*2
# with penalty='penalty' applied, solution is:
xs = [1.22797136, 4.24537337]
ys = -0.09582504

from mystic.symbolic import generate_constraint, generate_solvers, solve
from mystic.symbolic import generate_penalty, generate_conditions

equations = """
x0**2 - x1 + 1.0 <= 0.0
1.0 - x0 + (x1 - 4)**2 <= 0.0
"""
#cf = generate_constraint(generate_solvers(solve(equations))) #XXX: inequalities
pf = generate_penalty(generate_conditions(equations), k=1e12)

from mystic.constraints import as_constraint

cf = as_constraint(pf)



if __name__ == '__main__':

    from mystic.solvers import buckshot
    from mystic.math import almostEqual

    result = buckshot(objective, 2, 40, bounds=bounds, penalty=pf, disp=False, full_output=True)

    assert almostEqual(result[0], xs, tol=1e-2)

               -x[0] + 2*x[1] <= 4

    where:  0 <= x[0] <= inf
            0 <= x[1] <= 4
"""
import numpy as np
import mystic.symbolic as ms
import mystic.solvers as my
import mystic.math as mm

# generate constraints and penalty for a linear system of equations 
A = np.array([[1, 1],[-1, 2]])
b = np.array([7,4])
eqns = ms.linear_symbolic(G=A, h=b)
cons = ms.generate_constraint(ms.generate_solvers(ms.simplify(eqns)))
pens = ms.generate_penalty(ms.generate_conditions(eqns), k=1e3)
bounds = [(0., None), (0., 4.)]

# get the objective
def objective(x):
  x = np.asarray(x)
  return x[0]**2 + 4*x[1]**2 - 32*x[1] + 64

x0 = np.random.rand(2)

# compare against the exact minimum
xs = np.array([2., 3.])
ys = objective(xs)


sol = my.fmin_powell(objective, x0, constraint=cons, penalty=pens, disp=False,

def pf(len=3):
    return generate_penalty(generate_conditions(equations(len)))

T + H + E + M + E - 72 == 0
V + I + O + L + I + N - 100 == 0
W + A + L + T + Z - 34 == 0
"""
var = list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')
#NOTE: FOR A MORE DIFFICULT PROBLEM, COMMENT OUT THE FOLLOWING 5 LINES
bounds[0] = (5,5)    # A
bounds[4] = (20,20)  # E
bounds[8] = (25,25)  # I
bounds[14] = (10,10) # O
bounds[20] = (1,1)   # U

from mystic.constraints import unique, near_integers, has_unique

from mystic.symbolic import generate_penalty, generate_conditions
pf = generate_penalty(generate_conditions(equations,var),k=1)
from mystic.constraints import as_constraint
cf = as_constraint(pf)
from mystic.penalty import quadratic_equality

@quadratic_equality(near_integers)
@quadratic_equality(has_unique)
def penalty(x):
    return pf(x)

from numpy import round, hstack, clip
def constraint(x):
    x = round(x).astype(int) # force round and convert type to int
    x = clip(x, 1,nletters)  #XXX: hack to impose bounds
    x = unique(x, list(range(1,nletters+1)))
    return x