06/29/17

# Solving the Fish Riddle with JuMP

Recently I came across a nice Ted-Ed video presenting a Fish Riddle.

I thought it would be fun to try solving it using Julia’s award winning JuMP package. Before we get started, please watch the above video-you might want to pause at 2:24 if you want to solve it yourself.

To attempt this problem in Julia, you will have to install the JuMP package.

julia> Pkg.add("JuMP")

JuMP provides an algebraic modeling language for dealing with mathematical optimization problems. Basically, that allows you to focus on describing your problem in a simple syntax and it would then take care of transforming that description in a form that can be handled by any number of solvers. Those solvers can deal with several types of optimization problems, and some solvers are more generic than others. It is important to pick the right solver for the problem that you are attempting.

The problem premises are:
1. There are 50 creatures in total. That includes sharks outside the tanks and fish
2. Each SECTOR has anywhere from 1 to 7 sharks, with no two sectors having the same number of sharks.
3. Each tank has an equal number of fish
4. In total, there are 13 or fewer tanks
5. SECTOR ALPHA has 2 sharks and 4 tanks
6. SECTOR BETA has 4 sharsk and 2 tanks
We want to find the number of tanks in sector GAMMA!

Here we identify the problem as mixed integer non-linear program (MINLP). We know that because the problem involves an integer number of fish tanks, sharks, and number of fish inside each tank. It also non-linear (quadratic to be exact) because it involves multiplying two two of the problem variables to get the total number or creatures. Looking at the table of solvers in the JuMP manual. pick the Bonmin solver from AmplNLWriter package. This is an open source solver, so installation should be hassle free.

julia> Pkg.add("AmplNLWriter")

We are now ready to write some code.

using JuMP, AmplNLWriter

# Solve model
m = Model(solver=BonminNLSolver())

# Number of fish in each tank
@variable(m, n>=1, Int)

# Number of sharks in each sector
@variable(m, s[i=1:3], Int)

# Number of tanks in each sector
@variable(m, nt[i=1:3]>=0, Int)

@constraints m begin
# Constraint 2
sharks[i=1:3], 1 <= s[i] <= 7
numfish[i=1:3], 1 <= nt[i]
# Missing uniqueness in restriction
# Constraint 4
sum(nt) <= 13
# Constraint 5
s[1] == 2
nt[1] == 4
# Constraint 6
s[2] == 4
nt[2] == 2
end

# Constraints 1 & 3
@NLconstraint(m, s[1]+s[2]+s[3]+n*(nt[1]+nt[2]+nt[3]) == 50)

# Solve it
status = solve(m)

sharks_in_each_sector=getvalue(s)
fish_in_each_tank=getvalue(n)
tanks_in_each_sector=getvalue(nt)

@printf("We have %d fishes in each tank.\n", fish_in_each_tank)
@printf("We have %d tanks in sector Gamma.\n",tanks_in_each_sector[3])
@printf("We have %d sharks in sector Gamma.\n",sharks_in_each_sector[3])

In that representation we could not capture the restriction that “no two sectors having the same number of sharks”. We end up with the following output:

We have 4 fishes in each tank.
We have 4 tanks in sector Gamma.
We have 4 sharks in sector Gamma.

Since the problem domain is limited, we can possible fix that by adding a constrain that force the number of sharks in sector Gamma to be greater than 4.

@constraint(m,s[3]>=5)

This will result in an answer that that does not violate any of the stated constraints.

We have 3 fishes in each tank.
We have 7 tanks in sector Gamma.
We have 5 sharks in sector Gamma.

However, this seems like a bit of kludge. The proper way go about it is represent the number of sharks in the each sector as binary array, with only one value set to 1.

# Number of sharks in each sector
@variable(m, s[i=1:3,j=1:7], Bin)

We will have to modify our constraint block accordingly

@constraints m begin
# Constraint 2
sharks[i=1:3], sum(s[i,:]) == 1
u_sharks[j=1:7], sum(s[:,j]) <=1 # uniquness
# Constraint 4
sum(nt) <= 13
# Constraint 5
s[1,2] == 1
nt[1] == 4
# Constraint 6
s[2,4] == 1
nt[2] == 2
end

We invent a new variable array st to capture the number of sharks in each sector. This simply obtained by multiplying the binary array by the vector $[1,2,\ldots,7]^\top$

@variable(m,st[i=1:3],Int)
@constraint(m, st.==s*collect(1:7))

We rewrite our last constraint as

# Constraints 1 & 3
@NLconstraint(m, st[1]+st[2]+st[3]+n*(nt[1]+nt[2]+nt[3]) == 50)

After the model has been solved, we extract our output for the number of sharks.

sharks_in_each_sector=getvalue(st)

…and we get the correct output.

This problem might have been an overkill for using a full blown mixed integer non-linear optimizer. It can be solved by a simple table as shown in the video. However, we might not alway find ourselves in such a fortunate position. We could have also use mixed integer quadratic programming solver such as Gurobi which would be more efficient for that sort of problem. Given the small problem size, efficiency hardly matters here.

06/24/17

# The blinded eye of the mind

We humans have always seemed to have sort of ineptitude of not being able to notice the great diversity of life on our planet. However, any feeling of emptiness is merely an illusion to a blind mind. Through the course of human evolution, our minds and our vast intellect  were the most essential tools for the survival of our species.

But if you were to condense Earth’s entire existence into 24hrs, then we would have only been on this planet for a mere… 3 seconds !

Our minds have enabled us to think deeply about our purpose in life and to do incredible things. But we have become so blinded by the power we have that we have not looked  at what insidious and irrevocable damage we have caused to the balance nature.

Do you know what is the most dangerous flaw our human minds have kept for the last 200,000 of our existence?…

Not knowing where to look or whom to blame for our mistakes. We have just become so ignorant, so blinded, and so utterly self-centered about only our well-being that we have not noticed issues like increasing famine cases around Nigeria and  Somalia, or the rapid occurrences of some of the largest storms ever seen.  So pernicious we have become to ourselves. How have we not noticed this catastrophic conundrum at all ?

Sure we have accomplished some of the most incredible feats of innovation in the natural world but we have not at all been wise enough to contemplate on the consequences our triumph cost. It’s good to be smart but not too smart for you own good.

We have been under the highly foolish illusion that what our minds and ideas have enabled us to do  would not exponentially changing many highly essential factors for the greater good of our species and all life on Earth.

We have believed for so long that supposedly “small” problems that don’t fall into our best interests at heart should not be  look at seriously for even the smallest or seemingly “simple” could bring the down fall of us all.

Do you want what is the deadliest disease ever?

Tuberculosis?

Measles?

Ebola ?

AIDs?  *It’s a single, influential idea which makes us blind to what heinous and atrocious things happen by every passing day.

Is it highly ethical of us as a species to fall into the hands of such a vicious free-thinking paralytic? Would it help to just think that with a flick of a finger all of the problems that this crisis that we have fallen into would be solved ?

An idea can be very deadly by it self  but there is a much more vicious motivator for simply ignoring the worsening problems we have caused.

We are silenced by payment of money. A bribery which many us might  lightly take. Our minds can be easily persuaded into doing any thing  by the motivation of payment or wealth.  If we are given everything we could possibly need it’s still not enough to content us.

Our greedy and cunning simply want more of this pleasurable wealth. Some people may see a hundred pound not in there hands from there work.  I think we have just become slaves to a printed peace of paper.

Such as a humming bird and the nectar which

06/19/17

This is my first post on the blog. I love the sea.  This picture was taken when I was on vacation in the Red Sea.

06/12/17

# Reading DataFrames with non-UTF8 encoding in Julia

Recently I ran into problem where I was trying to read a CSV files from a Scandinavian friend into a DataFrame. I was getting errors it could not properly parse the latin1 encoded names.

I tried running

using DataFrames
dataT=readtable("example.csv", encoding=:latin1)

but the got this error

ArgumentError: Argument 'encoding' only supports ':utf8' currently.

The solution make use of (StringEncodings.jl)[https://github.com/nalimilan/StringEncodings.jl] to wrap the file data stream before presenting it to the readtable function.

f=open("example.csv","r")
s=StringDecoder(f,"LATIN1", "UTF-8")
close(s)
close(f)

The StringDecoder generates an IO stream that appears to be utf8 for the readtable function.

06/12/17

# Tupper’s self-referential formula in Julia

I was surprised when I came across on Tupper’s formula on twitter. I felt the compulsion to implement it in Julia.

The formula is expressed as

${1\over 2} < \left\lfloor \mathrm{mod}\left(\left\lfloor {y \over 17} \right\rfloor 2^{-17 \lfloor x \rfloor - \mathrm{mod}(\lfloor y\rfloor, 17)},2\right)\right\rfloor$

and yields bitmap facsimile of itself.

In [1]:
k=big"960 939 379 918 958 884 971 672 962 127 852 754 715 004 339 660 129 306 651 505 519 271 702 802 395 266 424 689 642 842 174 350 718 121 267 153 782 770 623 355 993 237 280 874 144 307 891 325 963 941 337 723 487 857 735 749 823 926 629 715 517 173 716 995 165 232 890 538 221 612 403 238 855 866 184 013 235 585 136 048 828 693 337 902 491 454 229 288 667 081 096 184 496 091 705 183 454 067 827 731 551 705 405 381 627 380 967 602 565 625 016 981 482 083 418 783 163 849 115 590 225 610 003 652 351 370 343 874 461 848 378 737 238 198 224 849 863 465 033 159 410 054 974 700 593 138 339 226 497 249 461 751 545 728 366 702 369 745 461 014 655 997 933 798 537 483 143 786 841 806 593 422 227 898 388 722 980 000 748 404 719"
setprecision(BigFloat,10000);


In the above, the big integer is the magic number that lets us generate the image of the formula. I also need to setprecision of BigFloat to be very high, as rounding errors using the default precision does not get us the desired results. The implementation was inspired by the one in Python, but I see Julia a great deal more concise and clearer.

In [2]:
function tupper_field(k)
field=Array{Bool}(17,106)
for (ix,x) in enumerate(0.0:1:105.0), (iy,y) in enumerate(k:k+16)
field[iy,107-ix]=1/2<floor(mod(floor(y/17)*2^(-17*floor(x)-mod(floor(y),17)),2))
end
field
end

In [3]:
f=tupper_field(k);
using Images
img = colorview(Gray,.!f)

Out[3]:

I just inverted the boolean array here to get the desired bitmap output.