# Solving the Fish Riddle with JuMP

Re-posted from: http://perfectionatic.org/?p=438

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.

# Exploring Fibonacci Fractions with Julia

Re-posted from: http://perfectionatic.org/?p=367

Recently, I came across a fascinating blog and video from Mind you Decisions. It is about how a fraction
$\frac{1}{999,999,999,999,999,999,999,998,,999,999,999,999,999,999,999,999}$
would show the Fibonacci numbers in order when looking at its decimal output.

On a spreadsheet and most standard programming languages, such output can not be attained due to the limited precision for floating point numbers. If you try this on R or Python, you will get an output of 1e-48.
Wolfram alpha,however,allows arbitrary precision.

In Julia by default we get a little better that R and Python

julia> 1/999999999999999999999998999999999999999999999999
1.000000000000000000000001000000000000000000000002000000000000000000000003000002e-48

julia> typeof(ans)
BigFloat

We observe here that we are getting the first few Fibonacci numbers $1, 1, 2, 3$. We need more precision to get more numbers. Julia has arbitrary precision arithmetic baked into the language. We can crank up the precision of the BigFloat type on demand. Of course, the higher the precision, the slower the computation and the greater the memory we use. We do that by setprecision.

julia> setprecision(BigFloat,10000)
10000

Reevaluating, we get

julia> 1/999999999999999999999998999999999999999999999999
1.00000000000000000000000100000000000000000000000200000000000000000000000300000000000000000000000500000000000000000000000800000000000000000000001300000000000000000000002100000000000000000000003400000000000000000000005500000000000000000000008900000000000000000000014400000000000000000000023300000000000000000000037700000000000000000000061000000000000000000000098700000000000000000000159700000000000000000000258400000000000000000000418100000000000000000000676500000000000000000001094600000000000000000001771100000000000000000002865700000000000000000004636800000000000000000007502500000000000000000012139300000000000000000019641800000000000000000031781100000000000000000051422900000000000000000083204000000000000000000134626900000000000000000217830900000000000000000352457800000000000000000570288700000000000000000922746500000000000000001493035200000000000000002415781700000000000000003908816900000000000000006324598600000000000000010233415500000000000000016558014100000000000000026791429600000000000000043349443700000000000000070140873300000000000000113490317000000000000000183631190300000000000000297121507300000000000000480752697600000000000000777874204900000000000001258626902500000000000002036501107400000000000003295128009900000000000005331629117300000000000008626757127200000000000013958386244500000000000022585143371700000000000036543529616200000000000059128672987900000000000095672202604100000000000154800875592000000000000250473078196100000000000405273953788100000000000655747031984200000000001061020985772300000000001716768017756500000000002777789003528800000000004494557021285300000000007272346024814100000000011766903046099400000000019039249070913500000000030806152117012900000000049845401187926400000000080651553304939300000000130496954492865700000000211148507797805000000000341645462290670700000000552793970088475700000000894439432379146400000001447233402467622100000002341672834846768500000003788906237314390600000006130579072161159100000009919485309475549700000016050064381636708800000025969549691112258500000042019614072748967300000067989163763861225800000110008777836610193100000177997941600471418900000288006719437081612000000466004661037553030900000754011380474634642900001220016041512187673800001974027421986822316700003194043463499009990500005168070885485832307200008362114348984842297700013530185234470674604900021892299583455516902600035422484817926191507500057314784401381708410100092737269219307899917600150052053620689608327700242789322839997508245300392841376460687116573000635630699300684624818301028472075761371741391301664102775062056366209602692574850823428107600904356677625885484473810507049252476708912581411411405930102594397055221918455182579303309636633329861112681897706691855248316295261201016328488578177407943098723020343826493703204299739348832404671111147398462369176231164814351698201718008635835925499096664087184867000739850794865805193502836665349891529892378369837405200686395697571872674070550577925589950242511475751264321287522115185546301842246877472357697022053e-48

That is looking much better. However it we be nice if we could extract the Fibonacci numbers that are buried in that long decimal. Using the approach in the original blog. We define a function

y(x)=one(x)-x-x^2

and calculate the decimal

a=1/y(big"1e-24")

Here we use the non-standard string literal big"..." to insure proper interpretation of our input. Using BigFloat(1e-24)) would first construct at floating point with limited precision and then do the conversion. The initial loss of precision will not be recovered in the conversion, and hence the use of big. Now we extract our Fibonacci numbers by this function

function extract_fib(a)
x=string(a)
l=2
fi=BigInt[]
push!(fi,1)
for i=1:div(length(x)-24,24)
j=parse(BigInt,x[l+1+(i-1)*24:l+i*24])
push!(fi,j)
end
fi
end

Here we first convert our very long decimal number of a string and they we exploit the fact the Fibonacci numbers occur in blocks that 24 digits in length. We get out output in an array of BigInt. We want to compare the output with exact Fibonacci numbers, we just do a quick and non-recursive implementation.

function fib(n)
f=Vector{typeof(n)}(n+1)
f[1]=f[2]=1;
for i=3:n+1
f[i]=f[i-1]+f[i-2]
end
f
end

Now we compare…

fib_exact=fib(200);
fib_frac=extract_fib(a);
for i in eachindex(fib_frac)
println(fib_exact[i], " ", fib_exact[i]-fib_frac[i])
end

We get a long sequence, we just focused here on when the discrepancy happens.

...
184551825793033096366333 0
298611126818977066918552 0
483162952612010163284885 0
781774079430987230203437 -1
1264937032042997393488322 999999999999999999999998
2046711111473984623691759 1999999999999999999999997
...

The output shows that just before the extracted Fibonacci number exceeds 24 digits, a discrepancy occurs. I am not quite sure why, but this was a fun exploration. Julia allows me to do mathematical explorations that would take one or even two orders of magnitude of effort to do in any other language.

# Mathematical Sandpiles

Re-posted from: http://perfectionatic.org/?p=320

I can across an exotic number system known as a sandpile on Numberphile video.

Exploiting Julia fantastic type system and very cool index handling, I put implements that sandpile addition operator demonstrated in the Numberphile video.

The notebook can be found here.