Announcing support for complex-domain linear programs in Convex.jl

I am pleased to announce the support for complex-domain linear programs (LPs) in Convex.jl. As one of the Google Summer of Code students under The Julia Language, I had proposed to implement the support for complex semidefinite programming. In the first phase of project, I started by tackling the problem of complex-domain LPs where in first subphase, I had announced the support for complex coefficients during JuliaCon’16 and now I take this opportunity to announce the support for complex variables in LPs.

Complex-domain LPs consist of a real linear objective function, real linear inequality constraints, and real and complex linear equality constraints.

In order to enable complex-domain LPs, we came up with these ideas:

1. We redefined the conic_form! of every affine atom to accept complex arguments.
2. Every complex variable z was internally represented as z = z1 + i*z2, where z1 and z2 are real.
3. We introduced two new affine atoms real and imag which return the real and the imaginary parts of the complex variable respectively.
4. transpose and ctranspose perform differently on complex variables so a new atom CTransposeAtom was created.
5. A complex-equality constraint RHS = LHS can be decomposed into two corresponding real equalities constraint real(RHS) = real(LHS) and imag(RHS) = imag(LHS)

After above changes were made to the codebase, we wrote two use cases to demonstrate the usability and the correctness of our idea which I am presenting below:

# Importing Packages
Pkg.clone("https://github.com/Ayush-iitkgp/Convex.jl/tree/gsoc2")
using Convex

# Complex LP with real variable
n = 10 # variable dimension (parameter)
m = 5 # number of constraints (parameter)
xo = rand(n)
A = randn(m,n) + im*randn(m,n)
b = A * xo
# Declare a real variable
x = Variable(n)
p1 = minimize(sum(x), A*x == b, x>=0)
# Notice A*x==b is complex equality constraint
solve!(p1)
x1 = x.value

# Let's now solve by decomposing complex equality constraint into the corresponding real and imaginary part.
p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0)
solve!(p2)
x2 = x.value
x1==x2 # should return true

# Let's now consider an example using a complex variable
# Complex LP with complex variable
n = 10 # variable dimension (parameter)
m = 5 # number of constraints (parameter)
xo = rand(n)+im*rand(n)
A = randn(m,n) + im*randn(m,n)
b = A * xo

# Declare a complex variable
x = ComplexVariable(n)
p1 = minimize(real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0)
solve!(p1)
x1 = x.value

xr = Variable(n)
xi = Variable(n)
p2 = minimize(sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0)
solve!(p2)
x1== xr.value + im*xi.value # should return true


List of all the affine atoms are as follows:

1. addition, substraction, multiplication, division
2. indexing and slicing
3. k-th diagonal of a matrix
4. construct diagonal matrix
5. transpose and ctranspose
6. stacking
7. sum
8. trace
9. conv
10. real and imag

Now, I am working towards implementing complex-domain second order conic programming. Meanwhile, I invite the Julia community to play around with the complex-domain LPs. The link to the development branch is here.

Looking forward to your suggestions!

Special thanks to my mentors Madeleine Udell and Dvijotham Krishnamurthy!

Introducing DifferentialEquations.jl

Differential equations are ubiquitous throughout mathematics and the sciences. In fact, I myself have studied various forms of differential equations stemming from fields including biology, chemistry, economics, and climatology. What was interesting is that, although many different people are using differential equations for many different things, pretty much everyone wants the same thing: to quickly solve differential equations in their various forms, and make some pretty plots to describe what happened.

The goal of DifferentialEquations.jl is to do exactly that: to make it easy solve differential equations with the latest and greatest algorithms, and put out a pretty plot. The core idea behind DifferentialEquations.jl is that, while it is easy to describe a differential equation, they have such diverse behavior that experts have spent over a century compiling different ways to think about and handle differential equations. Most users will want to just brush past all of the talk about which algorithms simply ask: “I have this kind of differential equation. What does the solution look like?”

DifferentialEquations.jl’s User Interface

To answer that question, the user should just have to say what their problem is, tell the computer to solve it, and then tell the computer to plot it. In DifferentialEquations.jl, we use exactly those terms. Let’s look at an Ordinary Differential Equation (ODE): the linear ODE . It is described as the function

$y^\prime = \alpha y$

To use DifferentialEquations.jl, you first have to tell the computer what kind of problem you have, and what your data is for the problem. Recall the general ordinary differential equation is of the form

$y^\prime = f(y)$

and initial condition $u_0$, so in this case, we have an ODE with data $f(y)=\alpha y$ and $u_0$. DifferentialEquations.jl is designed as a software for a high-level language, Julia. There are many reasons for this choice, but the one large reason is its type system and multiple dispatch. For our example, we tell the machine what type of problem we have by building a DEProblem type. The code looks like this:

using DifferentialEquations
alpha = 0.5 #Setting alpha to 1/2
f(y,t) = alpha*y
u0 = 1.5
prob = ODEProblem(f,u0)

where prob contains everything about our problem. You can then tell the computer to solve it and give you a plot by, well, solve and plot:

timespan = [0,1] # Solve from time = 0 to time = 1
sol = solve(prob,timespan) # Solves the ODE
plot(sol) # Plots the solution using Plots.jl

And that’s the key idea: the user should simply have to tell the program what the problem is, and the program should handle the details. That doesn’t mean that the user won’t have access to to all of the details. For example, we can control the solver and plotters in more detail, using something like

sol = solve(prob,alg=:RK4) # Unrolled and optimzed RK4
plot(sol,lw=3) # All of the Plots.jl attributes are available

However, in many cases a user may approach the problem for which they don’t necessarily know anything about the algorithms involved in approximating the problem, and so obfuscating the API with these names is simply confusing. One place where this occurs is solving stochastic differential equations (SDEs). These have been recently growing in popularity in many of the sciences (especially systems biology) due to their added realism and their necessity when modeling rare and random phenomena. In DifferentialEquations.jl, you can get started by simply knowing that an SDE problem is defined by the functions $f$ and $g$ in the form

$dX_t = f(X_t,t)dt + g(X_t,t)dW_t,$

with initial condition $u_0$, and so the steps for defining and solving the linear SDE is

g(u,t) = 0.3u
prob = SDEProblem(f,g,u0)
sol = solve(prob,timespan)
plot(sol)

If you wish to dig into the manual, you will see that the default solver that is used is a Rossler-SRI type of method and will (soon) have adaptivity which is complex enough to be a numerical analysis and scientific computing research project. And you can dig into the manual to find out how to switch to different solvers, but the key idea is that you don’t have to. Everything is simply about defining a problem, and asking for solutions and plots.

And that’s it. For more about the API, take a look at the documentation or the tutorial IJulia notebooks. What I want to discuss is why I believe this is the right choice, where we are, and where we can go with it.

What exactly does that give us?

Julia was created to solve the many-language problem in scientific computing. Before people would have to write out the inner loops as C/Fortran, and bind it to a scripting language that was never designed with performance in mind. Julia has done extremely well as solving this problem via multiple-dispatch. Multiple dispatch is not just about ease of use, but it is also the core of what allows Julia to be fast . From a quote I am stealing from IRC: “Julia: come for the fast loops, stay for the type-system”.

In my view, the many-language problem always had an uglier cousin: the many-API problem. Every package has its own way of interacting with the user, and it becomes a burden to remember how all of them work. However, in Julia there seems to be a beautiful emergence of packages which solve the many-API problem via Julia’s multiple-dispatch and metaprogramming functionalities. Take for example Plots.jl. There are many different plotting packages in Julia. However, through Plots.jl, you can plot onto any “backend” (other plotting package) using just one API. You can mix and match plotting in PyPlot (matplotlib), GR, Plotly, and unicode. It’s all the same commands. Another example of this is JuMP. Its core idea is solver independence: you take your optimization problem, define the model in JuMP’s lingo, and then plug into many different solvers all by flipping a switch.

DifferentialEquations.jl is extending this idea to the realm of differential equations. By using the keyword alg=:ode45, the solver can call functions from ODE.jl. And changing it to alg=:dopri5, DifferentialEquations.jl will solve your ODE problem using the coveted dopri5 Fortran software. The complexity of learning and understanding many different APIs is no longer a requirement for trying different algorithms.

But why “Differential Equations”? Isn’t that broad?

Sure, there are packages for solving various types of differential equations, all specializing in one little part. But when I was beginning my PhD, quickly found that these packages were missing something. The different types of differential equations that we encounter are not different but many times embody the same problem: a PDE when discretized is a system of ODEs, the probability distribution of evolving SDEs is a PDE (a form of the Heat Equation), and all of the tools that they use to get good performance are the same. Indeed, many contemporary research questions can be boiled down to addressing the question: what happens if we change the type of differential equation? What happens if we add noise to our ODEs which describe population dispersal? What happens if we add to our model that RNA production is reading a delayed signal? Could we make this high-dimensional PDE computationally feasible via a Monte Carlo experiment combined with Feynman-Kac’s theorem?

Yet, our differential equations libraries are separate. Our PDEs are kept separate from our SDEs, while our delay equations hang out in their own world. Mixing and matching solvers requires learning complex APIs, usually large C/Fortran libraries with opaque function names. That is what DifferentialEquations.jl is looking to solve. I am building DifferentialEquations.jl as a hub for differential equations, the general sense of the term.

If you have defined an SDE problem, then via the Forward Kolmorogov equation there is a PDE associated to the SDE. In many cases like the Black-Scholes model, both the SDE and the PDE are canonical ways of looking at the same problem. The solver should translate between them, and the solver should handle both types of differential equations. With one API and the functionality for these contained within the same package, no longer are they separate entities to handle computationally.

Where are we currently?

DifferentialEquations.jl is still very young. Indeed, the project only started a few months ago, and during that time period I was taking 6 courses. However, the package already has a strong core, including

In fact, there are already a lot of features which are unique to DifferentialEquations.jl:

You may have been thinking, “but I am a numerical analyst. How could this program help me?”. DifferentialEquations.jl has a bunch of functionalities for quickly designing and testing algorithms. All of the DEProblems allow for one to give them the analytical solution, and the solvers will then automatically calculate the errors. Thus by using some simple macros, one can define new algorithms in just a few lines of code, test the convergence, benchmark times, and have the algorithm available as an alg option in no time (note: all of the ODE solvers were written in one morning!). Thus it is easy to define the loop, and the rest of the functionality will come by default. It’s both a great way to test algorithms, and share algorithms. Contributing will both help you and DifferentialEquations.jl!.

Where are we going?

I have big plans for DifferentialEquations.jl. For example:

• I will be rolling out an efficiency testing suite so that one could just specify the algorithms you’d like to solve a problem, and along what convergence axis (i.e. choose a few $\Delta t$s, or along changing tolerances), and it will output comparisons of the computational efficiencies and make some plots. It will be similar in functionality to the ConvergenceSimulation suite.
• Finite difference methods for Heat and Poisson equation. These are long overdue for the research I do.
• Changing the tableaus in ODEs and SDEs to StaticArrays so they are stack allocated. This has already been tested and works well on v0.5.
• Higher-order methods for parabolic SPDEs (a research project with promising results!).
• Blazing fast adaptivity for SDEs. (Once the paper I have submitted for it is published, it will be available. It’s already implemented!)
• High-stability high order methods for SDEs (another research project).
• Parallel methods. I have already implemented parallel (Xeon Phi) solvers and described them in previous blog posts. They simply need to be integrated into DifferentialEquations.jl. I would like to have native GPU solvers as well.
• Delay and algebraic differential equations.
• Wrapping more popular solvers. I’d like to add Sundials, LSODE, and PetsC to the list.
• A web interface via Escher.jl to define DEProblems and get the solution plots. I am looking to have this hosted as an XSEDE Gateway.

If you’d like to influence where this project is going, please file an issue on the Github repository. I am always open for suggestions.

I hope this gives you a good idea on what my plans are for DifferentialEquations.jl. Check out the documentation and give the package a whirl!

The post Introducing DifferentialEquations.jl appeared first on Stochastic Lifestyle.

I’m Moving to Gadfly

### Christina Lee

So far in this blog, all my plots have either used PyPlot, or the experimental GLVisualize. BUT, NO LONGER!

While writing up my research notes, I’ve been trying to improve my PyPlot graphs. I was hoping for them to look nicer; Julia’s PyPlot was saying no. We had a falling out over background colors and resolution. So, I have decided to move on to a plotting partner more willing to cope with my demands.

When I first started learning Julia, Gadfly seemed too unfamiliar, and control too tricky. I didn’t give it much further thought. After some digging, Gadfly has changed my mind, and I hope today I can change yours too. Right now, like Julia, Gadfly doesn’t have the maturity or documentation of PyPlot, but once it does, you better watch out.

Gadfly preferentially produces plots in SVG (scalable vector graphics) form, instead of the rasterized form of PNG or JPG favored by PyPlot. As displayed in the figure left, rasterized images are composed of colored squares, so don’t look the best when highly magnifified. SVG, which is based on XML, stores the position of vector elements and scales to high resolution without artifacts. Tools such as Inkscape also work in SVG, and can edit the individual pieces of a graph since they haven’t been rasterized together.

Let’s plot some graphs!

By Yug, modifications by 3247 (Unknown) [CC BY-SA 2.5 (http://creativecommons.org/licenses/by-sa/2.5)], via Wikimedia Commons
</div>

# Time to load our Packages
#Pkg.update()

using Compose

# Some test data
x=collect(1:.4:10);
y1=sin(x);
y2=cos(x);

# A simple plot
plot(x=x,y=y1)


Unlike PyPlot, we can give plot functions as well as points. After receiving a function, a max value, and a min value, Gadfly figures everything else out for itself.

For the function, we can pass in either an inbuilt function, a user defined function, or an anonymous function.

The brackets [ ] allow us to group multiple functions together in the plot.

# plotting function

function testfunction(x::Real)
return 1-x^2
end

plot([cos,testfunction,x->x^2],0,1)


But that’s just one set of points…
How do we plot more than one set of data? That’s where layers come in. Instead of using plot, we use layer, and set it as a variable. We can then plot those layers.

p1=layer(x=x,y=y1,Geom.point)
p2=layer(x=x,y=y2,Geom.point)
plot(p1,p2)


Different style of plotting: Geom

Now, what if we want to plot something other than lines?
Gadfly allows a wide variety of other ways to plot our data. We control this with the Geom (Geometry) object.

plot(x=x,y=y1,Geom.point,Geom.line)

plot(x=x,y=y1,ymin=y1-.1,ymax=y1+.1,Geom.point,Geom.errorbar)

plot(x=x,y=y1,ymin=y1-.1,ymax=y1+.1,Geom.line,Geom.ribbon)

plot(x=x,y=y1,Geom.bar)

plot(x=x,y=y1,Geom.step)

plot(x=y1,y=y2,Geom.path())


Edit Point and Line style

Disclaimer: I’ve learned how to do the next two things by reading the code, not the documentation.

It seems like they are newer and less developed features than everything else I’m discussing here. The syntax seems less polished than in other areas, so I believe it’s still under construction.

z=collect(0:.2:2)

xp=[z;z+.5;z+1;z+1.5;z+2;z+2.5;z+3;z+3.5]
yp=[z;z;z;z;z;z;z;z]
sh=[ones(z);2*ones(z);3*ones(z);4*ones(z);5*ones(z);6*ones(z);7*ones(z);8*ones(z)]

plot(x=xp,y=yp,shape=sh,Geom.point,Theme(default_point_size=5pt))

# or Compose.mm for smaller sizes
# These ratios and numbers changed around how ever you like

dash = 6 * Compose.cm
dot = 1.5 * Compose.cm
gap = 1 * Compose.cm

l1=layer(x=x,y=zeros(x),Geom.line,Theme(line_style=[dot]))
l2=layer(x=x,y=ones(x),Geom.line,Theme(line_style=[dash]))
l3=layer(x=x,y=2*ones(x),Geom.line,Theme(line_style=[gap]))
l4=layer(x=x,y=3*ones(x),Geom.line,Theme(line_style=[dot,dash,dot]))
l5=layer(x=x,y=4*ones(x),Geom.line,Theme(line_style=[gap,dash]))

plot(l1,l2,l3,l4,l5,Coord.Cartesian(ymin=-.5,ymax=4.5),Theme(grid_line_width=0pt))


Guide: Labeling Axes

Where Geom alters how we plot, Guide alters the labeling.
Guide ties in with Compose.jl through the Guide.annotate command, but that will take a tutorial in itself.

plot(x=x,y=y2,color=y2,
Guide.xlabel("xlabel"),Guide.ylabel("ylabel"),Guide.title("Something explanatory"),
Guide.colorkey("y2"))


Here’s something we can do with a combination of Guide and Scale. Using Guide, we can set where we specifically want our xticks to be, namely at multiples of $\pi$. But then, the labeling would write out some irrational number, making the plot look horrible. So, we create a function that takes in a number and outputs a string label for Scale to use.

function labelname(x::Real)
n=round(Int,x/π) #nearest integer*pi
if n==0
return "0"
elseif n==1
return "π"
end
return("$n π") end xmarks=[0,π,2π,3π] ymarks=[-1,0,1] plot(x=x,y=y2, Guide.yticks(ticks=ymarks), Guide.xticks(ticks=xmarks),Scale.x_continuous(labels=labelname))  Some other cool things we can do with Scale: * Automatically transform the axis according to log, log10,log2, asinh, or sqrt. * Write numbers in :plain, :scientific, :engineering, or :auto plot(x->exp(x),0,10,Scale.y_log,Scale.x_continuous(format=:scientific))  Themes I mostly chose Gadfly because of the control I could have with the Theme command. http://dcjones.github.io/Gadfly.jl/themes.html has a much more exhaustive list than what I will be demonstrating with here. # Solarized Colors that I like working with # http://ethanschoonover.com/solarized using Colors base03=parse(Colorant,"#002b36"); base02=parse(Colorant,"#073642"); base01=parse(Colorant,"#586e75"); base00=parse(Colorant,"#657b83"); base0=parse(Colorant,"#839496"); base1=parse(Colorant,"#839496"); base2=parse(Colorant,"#eee8d5"); base3=parse(Colorant,"#fdf6e3"); yellow=parse(Colorant,"#b58900"); orange=parse(Colorant,"#cb4b16"); red=parse(Colorant,"#dc322f"); magenta=parse(Colorant,"#d33682"); violet=parse(Colorant,"#6c71c4"); blue=parse(Colorant,"#268bd2"); cyan=parse(Colorant,"#3aa198"); green=parse(Colorant,"#859900");  plot(x=x,y=y1, Theme(highlight_width=0pt, # lose the white ring around the points default_point_size=3pt, # size of the dots default_color=magenta, # color of the dots background_color=base03, # ... background color ... grid_color=base2, # the lines minor_label_color=base2, # numbers on axes color key_label_color=base2, # color key numbers color key_title_color=cyan, # label of color key color major_label_color=cyan, # axes label and title color major_label_font_size=20pt, # font size for axes label and title minor_label_font_size=15pt, #font size for numbers on axes panel_opacity=1 #setting background to opaque );)  Coord: Setting the boundries The documentation states that you can change the range of the axes using Scale, but I’ve found it better to use this format to set my min and max values. plot(x=x,y=y1,Geom.line,Coord.Cartesian(ymin=0,ymax=1,xmin=0,xmax=10))  So far, we have covered seperately Guide, Themes, and partially Coord and Scale. Individually, each aspect doesn’t add too much clunkiness to the code. However, if we start to add everything together, then the plot function would look quite nasty. Luckily, just like layers for data points, we can put our modifiers into variables. Then we can simply call the variables in our plot function. This also helps for when we want to use one theme for a series of graphs. We can define the theme variable up top, and then include it in every graph there after, never having to declare it again. This helps me to keep my graphs uniform in style. function labelname(x::Real) n=round(Int,x/π) if n==0 return "0" elseif n==1 return "π" else return("$n π")
end
return("\$n π")
end

xticks=[0,π,2π,3π]
yticks=[-1,0,1]

data=layer(x=x,y=y1,ymin=y1-.1,ymax=y1+.1,Geom.line,Geom.ribbon)
f=layer(cos,0,3π)

yt=Guide.yticks(ticks=yticks)
xt=Guide.xticks(ticks=xticks)
xs=Scale.x_continuous(labels=labelname)

t= Theme(highlight_width=0pt,
default_point_size=3pt,
default_color=blue,
background_color=base03,
grid_color=base2,
minor_label_color=base2,
key_label_color=base2,
key_title_color=cyan,
major_label_color=cyan,
major_label_font_size=20pt,
minor_label_font_size=15pt,
panel_opacity=1)

xl=Guide.xlabel("time")
yl=Guide.ylabel("y axis")
GT=Guide.title("Important title")

plot(data,f,yt,xt,xs,t,xl,yl,GT)


Although we still have to call quite a few variables, this is a much simpler way of doing things.

Saving the Figure

Julia naturally saves to SVG (or SVGJS).
We have to specify the x and y dimensions of the plot, but since these images rescale so easily, we can choose some reasonable numbers.

p=plot(data,f,yt,xt,xs,t,xl,yl,GT)
draw(SVG("myplot.svg", 15cm, 9cm), p)


If you want to save to PNG, PDF, or PS, the package Cairo.jl provides that ability.

using Cairo
draw(PNG("myplot.png",15cm,9cm),p)


# Happy Plotting!

“julia

“