What is Objective coefficient for new variable - gurobi

In order to create new GRBVar, I need to provide Objective coefficient for new variable:
GRBVar var = model.addVar (double lowerBound,
double upperBound,
double obj, // objective coefficient
char type,
String name);
According to this example, the value can be set to 0. So I wonder what objective coefficient is.

Objective coefficient is the coefficient of the variable in your objective function. In the example you have given :
maximize x + y + 2 z
subject to x + 2 y + 3 z <= 4
x + y >= 1
x, y, z binary
your objective function is maximize x + y + 2 z
so Objective coefficients are
for x: 1
for y: 1
and for z: 2
While creating variables you can give coefficients arbitrary ( here they are as you said 0.0 )
// Create variables
GRBVar x = model.addVar(0.0, 1.0, 0.0, GRB.BINARY, "x");
GRBVar y = model.addVar(0.0, 1.0, 0.0, GRB.BINARY, "y");
GRBVar z = model.addVar(0.0, 1.0, 0.0, GRB.BINARY, "z");
But later you should set to the actual objective coefficients:
// Set objective: maximize x + y + 2 z
GRBLinExpr expr = new GRBLinExpr();
expr.addTerm(1.0, x);
expr.addTerm(1.0, y);
expr.addTerm(2.0, z);
model.setObjective(expr, GRB.MAXIMIZE);

Related

Formula to find intersection point of two lines

How to find an point where line1 and lin2 intersect, if both lines are defined by x,y,alpha where x,y are coordinates of a point on a line and alpha is the angle between the line and x=const?
I tried applying sine theorem but it yields two answers (triangles can be built on both sides of a line). I can check which point forms correct slope with one of the points, but that is ugly.
I can switch to y=ax+b representation, but then I have special cases that I have to worry about. Vertical and horizontal lines should be differently to avoid division by zero in 1/sin(alpha) and 1/cos(alpha) cases.
I am not looking for an implementation in a certain language, just a formula.
These questions are not relevant because they deal with finite line segments, NOT lines.
Given two points and two vectors, find point of intersection
How do you detect where two line segments intersect?
Suppose line 1 is defined by [x1, y1] and alpha1 and line 2 by [x2, y2] and alpha2.
Suppose k1 = tan(alpha1) and k2 = tan(alpha2).
Then the formula for the x-coordinate of the intersection is
x = (y2 - y1 + k1 * x1 - k2 * x2) / (k1 - k2)
Note: Function tan is undefined for angles pi / 2 + k * pi (where k is an arbitrary integer), so:
if k1 is undefined, then x = x1 and y = y2 + k2 * (x1 - x2)
if k2 is undefined, then x = x2 and y = y1 + k1 * (x2 - x1)
(both are practically the same with exchange of indices 1 <--> 2).
For a line equation Y = aX + b, you can calculate a = tan(alpha).
So if line1 is defined as x, y and alpha, the equation is Y = tan(alpha) * X + b.
Now to find b, you need a point on your line. This point has coordinates (x, y).
y = ax + b
b = y - ax
So you line equation is:
Y = tan(alpha) * X + (y - tan(alpha) * x)
Now you only have to solve the lines equation:
Y = a1 * X + b1
Y = a2 * X + b2
Which is:
a1 * X + b1 = a2 * X + b2
(a1 - a2) * X = b2 - b1
X = (b2 - b1) / (a1 - a2)
Right now you can calculate Y too.
So if we replace, we obtain:
X = ((y2 - tan(alpha2) * x2) - (y1 - tan(alpha1) * x1)) / (tan(alpha1) - tan(alpha2)
Simplified:
X = (y2 - y1 - tan(alpha2) * x2 + tan(alpha1) * x1)) / (tan(alpha1) - tan(alpha2)
And then:
Y = tan(alpha1) * X + (y - tan(alpha1) * x

Changing equation to extend integer overflow

Problem definition
I have a 2d scenario represented by a set of points (x,y), where x and y are 64bit integers. Both x and y values are in range of [0, R] (minimum possible x value is 0 and maximum value is R. Same rule applies for y).
At some point in my code i perform the following operation.
// Returns true if point P is inside circle defined by points A, B, C.
private bool InCircle(Long2 A, Long2 B, Long2 C, Long2 P)
{
Long2 AP = P - A;
Long2 BP = P - B;
Long2 CP = P - C;
return AP.SquareMagnitude * Long2.Det(BP, CP) + BP.SquareMagnitude * Long2.Det(CP, AP) + CP.SquareMagnitude * Long2.Det(AP, BP) > 0;
}
// Notes:
SquareMagnitude = x * x + y * y
Long2.Det(a, b) = a.x * b.y - a.y * b.x
I don't want this operation to fail, so the integers in it must not overflow. Maximum values occur when:
A = (0, 0)
B = (R, 0)
C = (0, R)
P = (R, R)
These points lead to:
AP = (R, R)
BP = (0, R)
CP = (R, 0)
If these values run on the method, we get:
AP.SquareMagnitude = R² + R²
Long2.Det(BP, CP) = R²
BP.SquareMagnitude = R²
Long2.Det(CP, AP) = R²
CP.SquareMagnitude = R²
Long2.Det(AP, BP) = R²
And the result becomes:
(2 * R² * R²) + (R² * R²) + (R² * R²) > 0
Which is the same as:
4 * (R^4) > 0
If we don't want a 64bit integer to overflow, then:
4 * (R^4) <= (2^63)
Thus, maximum R value is 2^15 = 32768. (It's actually 2^15.25, but we'll round it to 2^15).
Question
Is there any way to increase maximum R value without getting overflow using 64bit integers? Something like breaking the equation down into smaller values and comparing them separately.
My thoughts so far
I thought of changing the equation to:
return AP.SquareMagnitude * Long2.Det(BP, CP) > -BP.SquareMagnitude * Long2.Det(CP, AP) - CP.SquareMagnitude * Long2.Det(AP, BP);
Which would lead to:
2 * R² * R² > -(R² * R²) - (R² * R²)
2 * (R^4) <= 2^63
R <= 2^15.5 (which still gets rounded to 2^15. Not good enough)
Besides, i don't know if it's safe to do this way. Is it?
Using integers is intended and output accuracy is a must.
I'm aware of this library:
https://www.cs.cmu.edu/~quake/robust.html
But i want to use integers, not floats. The method being used here is robust within the integers overflow limitation. It is only desired to extend this overflow limitation.
So i thought of running checked integer arithmetic and if it overflows, convert the 64bit integers to decimals.
// Returns true if point P is inside circle defined by points A, B, C.
private bool InCircle(Long2 A, Long2 B, Long2 C, Long2 P)
{
Long2 AP = P - A;
Long2 BP = P - B;
Long2 CP = P - C;
try
{
long a = checked(AP.SquareMagnitude * Long2.Det(BP, CP));
long b = checked(BP.SquareMagnitude * Long2.Det(CP, AP));
long c = checked(CP.SquareMagnitude * Long2.Det(AP, BP));
long lhs = checked(checked(a + b) + c);
return lhs > 0;
}
catch (System.OverflowException)
{
decimal a1 = new decimal(AP.SquareMagnitude);
decimal a2 = new decimal(Long2.Det(BP, CP));
decimal b1 = new decimal(BP.SquareMagnitude);
decimal b2 = new decimal(Long2.Det(CP, AP);
decimal c1 = new decimal(CP.SquareMagnitude);
decimal c2 = new decimal(Long2.Det(AP, BP);
decimal lhs = a1 * a2 + b1 * b2 + c1 * c2;
return lhs > 0;
}
}
// Notes:
SquareMagnitude = x * x + y * y
Long2.Det(a, b) = a.x * b.y - a.y * b.x
Decimals can precisely represent integers < 2^96 if I'm not wrong. So to avoid precision loss and overflows:
4 * (R^4) <= 2^96
R <= 2^23 (It's actually 2^23.5)
And to avoid overflows in the integers operations, (R² + R²) is the largest operation:
R² + R² <= 2^63
2 * R² <= 2^63
R <= 2^31
Using this method, maximum R value becomes 2^23, which is great.
Another thought: would the code run faster on average if the method would not try to perform checked integer arithmetic and instead converting the values to decimals and going decimals math right from the start?
Is this going to work? Did I make something wrong?

maxima CAS - how to substitute variable for an expression?

In maxima, is there a way to apply variable substitutions for a subexpression? For example, replace instances of x+y with z.
subst works for the trivial case, but not for anything more than that.
(%i92) subst(x + y = foo, x + y);
(%o93) foo
(%i94) subst(x + y = foo, x + y + z);
(%o95) z + y + x
I think ratsubst has the effect you want.
(%i2) ratsubst(foo, x+y, x+y+z);
(%o2) z + foo

FORTRAN counter loop returns multiple iterations of the same value

First of all I am a complete novice to FORTRAN. With that said I am attempting to "build" a box, then randomly generate x, y, z coordinates for 100 atoms. From there, the goal is to calculate the distance between each atom, which becomes the value "r" of the Lennard-Jones potential energy equation. Then calculate the LJ potential, and finally sum the potential of the entire box. A previous question that I had asked about this project is here. The problem is that I get the same calculated value over and over and over again. My code is below.
program energytot
implicit none
integer, parameter :: n = 100
integer :: i, j, k, seed(12)
double precision :: sigma, r, epsilon, lx, ly, lz
double precision, dimension(n) :: x, y, z, cx, cy, cz
double precision, dimension(n*(n+1)/2) :: dx, dy, dz, LJx, LJy, LJz
sigma = 4.1
epsilon = 1.7
!Box length with respect to the axis
lx = 15
ly = 15
lz = 15
do i=1,12
seed(i)=i+3
end do
!generate n random numbers for x, y, z
call RANDOM_SEED(PUT = seed)
call random_number(x)
call random_number(y)
call random_number(z)
!convert random numbers into x, y, z coordinates
cx = ((2*x)-1)*(lx*0.5)
cy = ((2*y)-1)*(lx*0.5)
cz = ((2*z)-1)*(lz*0.5)
do j=1,n-1
do k=j+1,n
dx = ABS((cx(j) - cx(k)))
LJx = 4 * epsilon * ((sigma/dx(j))**12 - (sigma/dx(j))**6)
dy = ABS((cy(j) - cy(k)))
LJy = 4 * epsilon * ((sigma/dy(j))**12 - (sigma/dy(j))**6)
dz = ABS((cz(j) - cz(k)))
LJz = 4 * epsilon * ((sigma/dz(j))**12 - (sigma/dz(j))**6)
end do
end do
print*, dx
end program energytot
What exactly is your question? What do you want your code to do, and what does it do instead?
If you're having problems with the final print statement print*, dx, try this instead:
print *, 'dx = '
do i = 1, n * (n + 1) / 2
print *, dx(i)
end do
It seems that dx is too big to be printed without a loop.
Also, it looks like you're repeatedly assigning the array dx (and other arrays in the loop) to a single value. Try this instead:
i = 0
do j=1,n-1
do k=j+1,n
i = i + 1
dx(i) = ABS((cx(j) - cx(k)))
end do
end do
This way, each value cx(j) - cx(k) gets saved to a different element of dx, instead of overwriting previously saved values.
My new code goes something like this:
program energytot
implicit none
integer, parameter :: n = 6
integer :: i, j, k, seed(12)
double precision :: sigma, r, epsilon, lx, ly, lz, etot, pot, rx, ry, rz
double precision, dimension(n) :: x, y, z, cx, cy, cz
sigma = 4.1
epsilon = 1.7
etot=0
!Box length with respect to the axis
lx = 15
ly = 15
lz = 15
do i=1,12
seed(i)=i+90
end do
!generate n random numbers for x, y, z
call RANDOM_SEED(PUT = seed)
call random_number(x)
call random_number(y)
call random_number(z)
!convert random numbers into x, y, z coordinates
cx = ((2*x)-1)*(lx*0.5)
cy = ((2*y)-1)*(lx*0.5)
cz = ((2*z)-1)*(lz*0.5)
do j=1,n-1
do k=j+1,n
rx = (cx(j) - cx(k))
ry = (cy(j) - cy(k))
rz = (cz(j) - cz(k))
!Apply minimum image convention
rx=rx-lx*anint(rx/lx)
ry=ry-ly*anint(ry/ly)
rz=rz-lz*anint(rz/lz)
r=sqrt(rx**2+ry**2+rz**2)
pot=4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)
print*,pot
etot=etot+pot
end do
end do
print*, etot
end program energytot

Reflection of a Point over a Line

I have been looking at how to reflect a point in a line, and found this question which seems to do the trick, giving this formula to calculate the reflected point:
Given (x,y) and a line y = ax + c we want the point (x', y') reflected on the line.
Set d:= (x + (y - c)*a)/(1 + a^2)
Then x' = 2*d - x
and y' = 2*d*a - y + 2c
However there are two problems with this implementation for my needs:
My line is not described in the form y = ax + c (so I'd have to translate it, which is easy to do, but it means the process is slower).
What if a is infinity ie. a vertical line?
Is there a simple way to calculate (x', y'), the reflection of point (x, y) in a line, where the line is described by the two points (x1, y1) and (x2, y2)?
Edit:
I've found a formula which does this, but it seems as though it does not work with lines that look like they have equation y = x.
Here it is in actionscript:
public static function reflect(p:Point, l:Line):Point
{
// (l.sx, l.sy) = start of line
// (l.ex, l.ey) = end of line
var dx:Number = l.ex - l.sx;
var dy:Number = l.ey - l.sy;
if ((dx == 0) && (dy == 0))
{
return new Point(2 * l.sx - p.x, 2 * l.sy - p.y);
}
else
{
var t:Number = ((p.x - l.sx) * dx + (p.y - l.sy) * dy) / (dx * dx + dy * dy);
var x:Number = 2 * (l.sx + t * dx) - p.x;
var y:Number = 2 * (l.sy + t * dy) - p.y;
return new Point(x, y);
}
}
Does anyone have any idea where this formula goes wrong? I am still happy to take other solutions than the above formula - anything that'll work!
Find the normal vector from the point onto the line and add it twice to the point.
See wikipedia for the formula.
If you express everything in vectors you won't have the problems with an infinite slope.
;; line defined by ax+by+c=0
;; normal (a b) and distance from origin -c
(defun reflect-point-on-line (a b c px py)
(/ (+ (* a px)
(* b py)
c)
(sqrt (+ (expt a 2) (expt b 2)))))
#+nil ;; y-axis to (2 1)
(reflect-point-on-line 1 0 0 2 1) ;; => 2
#+nil ;; x-axis to (4 5)
(reflect-point-on-line 0 1 0 4 5) ;; => 5
If your line is 45 degrees... y = x...then your point (x1,y1) is one corner of a square, your line has two points on it that correspond to your point (x1 + distance to line, y1) and (x1, y1 plus distance to line) and the point you're looking for is the opposite corner, no?
If you find your distance, you should be able to add them to your coords and get your new point.

Resources