Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Appearance settings

Commit ba5b154

Browse filesBrowse files
yueminglyuemingl
committed
add paper examples
1 parent bfb1963 commit ba5b154
Copy full SHA for ba5b154

File tree

Expand file treeCollapse file tree

7 files changed

+187
-72
lines changed
Open diff view settings
Filter options
Expand file treeCollapse file tree

7 files changed

+187
-72
lines changed
Open diff view settings
Collapse file

‎src/symjava/domains/Domain.java‎

Copy file name to clipboardExpand all lines: src/symjava/domains/Domain.java
+10-7Lines changed: 10 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,9 @@ public static class CoordVarInfo {
2323
Double stepSize;
2424
Expr minBound;
2525
Expr maxBound;
26+
public String toString() {
27+
return minBound+" "+maxBound+" "+stepSize;
28+
}
2629
}
2730

2831
protected Map<Expr, CoordVarInfo> infoMap = new HashMap<Expr, CoordVarInfo>();
@@ -181,13 +184,13 @@ public Domain setConstraint(Expr logicalExpr) {
181184
return this;
182185
}
183186

184-
public Domain setConstraint(boolean dummy) {
185-
if(Ge.stackTop != null) {
186-
this.constraint = Ge.stackTop;
187-
Ge.stackTop = null;
188-
}
189-
return this;
190-
}
187+
// public Domain setConstraint(boolean dummy) {
188+
// if(Ge.stackTop != null) {
189+
// this.constraint = Ge.stackTop;
190+
// Ge.stackTop = null;
191+
// }
192+
// return this;
193+
// }
191194

192195
public Expr getConstraint() {
193196
return this.constraint;
Collapse file

‎src/symjava/examples/BlackScholes.java‎

Copy file name to clipboardExpand all lines: src/symjava/examples/BlackScholes.java
+38-3Lines changed: 38 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,8 +12,9 @@
1212
public class BlackScholes {
1313

1414
public static void main(String[] args) {
15-
test1(); //Example from QuantLib
16-
test2(); //Example from UCLA Statistics C183/C283: Statistical Models in Finance
15+
//test1(); //Example from QuantLib
16+
//test2(); //Example from UCLA Statistics C183/C283: Statistical Models in Finance
17+
peper_example();
1718
}
1819

1920

@@ -58,7 +59,7 @@ public static void test1() {
5859
BytecodeFunc blackScholesPrice = JIT.compile(new Expr[]{spot, strike, rd, rf, vol, tau, phi}, res);
5960
double price = blackScholesPrice.apply(100.0, 110.0, 0.002, 0.01, 0.1423, 0.5, 1);
6061

61-
System.out.println("Use Newtom method to recover the volatility by given the market data:");
62+
System.out.println("Use Newtom method to recover the volatility by giving the market data:");
6263
Expr[] freeVars = {vol};
6364
Expr[] params = {spot, strike, rd, rf, tau, phi}; //Specify the params in the given order
6465
Eq[] eq = new Eq[] {
@@ -144,4 +145,38 @@ public static void test2() {
144145
double[] constParams = new double[] {21, 20, 0.1, 0.25, 1.875};
145146
Newton.solve(eq, guess, constParams, 100, 1e-5);
146147
}
148+
149+
150+
public static void peper_example() {
151+
// Define symbols to construct the Black-Scholes formula
152+
Symbol spot = new Symbol("spot"); //spot price
153+
Symbol strike = new Symbol("strike"); //strike price
154+
Symbol rd = new Symbol("rd");
155+
Symbol rf = new Symbol("rf");
156+
Symbol vol = new Symbol("\\sigma"); //volatility
157+
Symbol tau = new Symbol("\\tau");
158+
Symbol phi = new Symbol("\\phi");
159+
160+
Expr domDf = exp(-rd*tau);
161+
Expr forDf = exp(-rf*tau);
162+
Expr fwd=spot*forDf/domDf;
163+
Expr stdDev=vol*sqrt(tau);
164+
//We use -10 instead of -oo for numerical computation
165+
double step = 1e-3;
166+
Domain I1 = Interval.apply(-10, phi*(log(fwd/strike)+0.5*pow(stdDev,2))/stdDev, z)
167+
.setStepSize(step);
168+
Domain I2 = Interval.apply(-10, phi*(log(fwd/strike)-0.5*pow(stdDev,2))/stdDev, z)
169+
.setStepSize(step);
170+
Expr cdf1 = Integrate.apply(exp(-0.5*pow(z,2)), I1)/sqrt(PI2);
171+
Expr cdf2 = Integrate.apply(exp(-0.5*pow(z,2)), I2)/sqrt(PI2);
172+
Expr res = phi*domDf*(fwd*cdf1-strike*cdf2);
173+
//Use Newtom method to recover the volatility by giving the market data
174+
Expr[] freeVars = {vol};
175+
Expr[] params = {spot, strike, rd, rf, tau, phi};
176+
Eq[] eq = new Eq[] { new Eq(res-0.1423, C0, freeVars, params) };
177+
// Use Newton's method to find the root
178+
double[] guess = new double[]{ 0.10 };
179+
double[] constParams = new double[] {100.0, 110.0, 0.002, 0.01, 0.5, 1};
180+
Newton.solve(eq, guess, constParams, 100, 1e-5);
181+
}
147182
}
Collapse file

‎src/symjava/examples/Example6.java‎

Copy file name to clipboardExpand all lines: src/symjava/examples/Example6.java
+45-29Lines changed: 45 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -31,37 +31,53 @@
3131
*/
3232
public class Example6 {
3333
public static void main(String[] args) {
34-
Func u = new Func("u", x, y);
35-
Func v = new Func("v", x, y);
36-
37-
// //Our PDE equation
38-
// Eq pde = new Eq(0.5*dot(grad(u), grad(v)) + 0.1*u*v, (x*x+y*y)*v);
39-
// //Read the mesh
40-
// Mesh2D mesh = new Mesh2D("mesh1", x, y);
41-
// mesh.readTriangleMesh("double_hex3.1.node", "double_hex3.1.ele");
42-
// solve(pde, mesh, null, "double_hex3.1.dat");
43-
44-
//Another PDE equation with Dirichlet condition
45-
WeakForm wf = new WeakForm(dot(grad(u), grad(v)), (-2*(x*x+y*y)+36)*v, u, v);
46-
//Eq pde2 = new Eq(u*v, (-2*(x*x+y*y)+36)*v);
47-
Mesh2D mesh2 = new Mesh2D("mesh2");
48-
//mesh2.readGridGenMesh("patch_triangle.grd");
49-
mesh2.readGridGenMesh("triangle.grd");
50-
//Mark boundary nodes
51-
double eps = 0.01;
52-
for(Node n : mesh2.nodes) {
53-
//if(1-Math.abs(n.coords[0])<eps || 1-Math.abs(n.coords[1])<eps || Math.abs(n.coords[0])<eps || Math.abs(n.coords[1])<eps )
54-
if(Math.abs(3-Math.abs(n.coords[0]))<eps || Math.abs(3-Math.abs(n.coords[1]))<eps)
55-
n.setType(1);
56-
}
57-
Map<Integer, Double> diri = new HashMap<Integer, Double>();
58-
diri.put(1, 0.0);
59-
//solve(pde2, mesh2, diri, "patch_triangle.dat");
60-
solve(wf, mesh2, diri, "triangle.dat");
61-
solve2(mesh2, diri, "triangle_hardcode.dat");
62-
34+
// Func u = new Func("u", x, y);
35+
// Func v = new Func("v", x, y);
36+
//
37+
//// //Our PDE equation
38+
//// Eq pde = new Eq(0.5*dot(grad(u), grad(v)) + 0.1*u*v, (x*x+y*y)*v);
39+
//// //Read the mesh
40+
//// Mesh2D mesh = new Mesh2D("mesh1", x, y);
41+
//// mesh.readTriangleMesh("double_hex3.1.node", "double_hex3.1.ele");
42+
//// solve(pde, mesh, null, "double_hex3.1.dat");
43+
//
44+
// //Another PDE equation with Dirichlet condition
45+
// WeakForm wf = new WeakForm(dot(grad(u), grad(v)), (-2*(x*x+y*y)+36)*v, u, v);
46+
// //Eq pde2 = new Eq(u*v, (-2*(x*x+y*y)+36)*v);
47+
// Mesh2D mesh2 = new Mesh2D("mesh2");
48+
// //mesh2.readGridGenMesh("patch_triangle.grd");
49+
// mesh2.readGridGenMesh("triangle.grd");
50+
// //Mark boundary nodes
51+
// double eps = 0.01;
52+
// for(Node n : mesh2.nodes) {
53+
// //if(1-Math.abs(n.coords[0])<eps || 1-Math.abs(n.coords[1])<eps || Math.abs(n.coords[0])<eps || Math.abs(n.coords[1])<eps )
54+
// if(Math.abs(3-Math.abs(n.coords[0]))<eps || Math.abs(3-Math.abs(n.coords[1]))<eps)
55+
// n.setType(1);
56+
// }
57+
// Map<Integer, Double> diri = new HashMap<Integer, Double>();
58+
// diri.put(1, 0.0);
59+
// //solve(pde2, mesh2, diri, "patch_triangle.dat");
60+
// solve(wf, mesh2, diri, "triangle.dat");
61+
// solve2(mesh2, diri, "triangle_hardcode.dat");
62+
peper_example();
6363
}
6464

65+
public static void peper_example() {
66+
Func u = new Func("u", x, y);
67+
Func v = new Func("v", x, y);
68+
WeakForm wf = new WeakForm(dot(grad(u), grad(v)), (-2*(x*x+y*y)+36)*v, u, v);
69+
Mesh2D mesh = new Mesh2D("Square");
70+
mesh.readGridGenMesh("triangle.grd");
71+
//Mark boundary nodes
72+
double eps = 0.01;
73+
for(Node n : mesh.nodes) {
74+
if(Math.abs(3-Math.abs(n.coords[0]))<eps || Math.abs(3-Math.abs(n.coords[1]))<eps)
75+
n.setType(1);
76+
}
77+
Map<Integer, Double> diri = new HashMap<Integer, Double>();
78+
diri.put(1, 0.0);
79+
solve(wf, mesh, diri, "triangle.dat");
80+
}
6581
public static void solve(WeakForm wf, Mesh2D mesh, Map<Integer, Double> dirichlet, String output) {
6682
System.out.println(String.format("PDE Weak Form: %s = %s", wf.lhs(), wf.rhs()));
6783

Collapse file

‎src/symjava/examples/NumericalIntegration.java‎

Copy file name to clipboardExpand all lines: src/symjava/examples/NumericalIntegration.java
+73-17Lines changed: 73 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
package symjava.examples;
22

3+
import symjava.relational.Ge;
34
import symjava.relational.Le;
45
import symjava.symbolic.*;
56
import static symjava.math.SymMath.*;
67
import static symjava.symbolic.Symbol.*;
78
import symjava.bytecode.BytecodeFunc;
89
import symjava.domains.Domain;
9-
import symjava.domains.Domain1D;
1010
import symjava.domains.Domain2D;
1111
import symjava.domains.Domain3D;
1212
import symjava.domains.DomainND;
@@ -17,20 +17,16 @@
1717
public class NumericalIntegration {
1818

1919
public static void main(String[] args) {
20-
//test_1D();
21-
//test_2D();
22-
//test_ND();
20+
// test_1D();
21+
// test_2D();
22+
// test_ND();
2323

2424
//Expr i = Integrate.apply(exp(pow(x,2)), Interval.apply(a, b).setStepSize(0.001));
2525
//BytecodeFunc fi = JIT.compile(new Expr[]{a,b}, i);
2626
//System.out.println(fi.apply(1,2));
2727

28-
Domain2D d = new Domain2D("D");
29-
d.setConstraint(
30-
x <=y & x <=z
31-
);
32-
System.out.println(d.getConstraint());
33-
28+
test_paper_example1();
29+
test_paper_example2();
3430
}
3531

3632
public static void test_1D() {
@@ -90,16 +86,16 @@ public static void test_2D() {
9086
.setBound(x, 0, 1)
9187
.setBound(y, 0, 1);
9288

93-
Expr ii = Integrate.apply(sin(sqrt(log(x+y+1))), omega);
94-
System.out.println(ii);
95-
BytecodeFunc f = JIT.compile(ii);
96-
System.out.println(f.apply());
89+
Expr I = Integrate.apply(sin(sqrt(log(x+y+1))), omega);
90+
System.out.println(I);
91+
BytecodeFunc fI = JIT.compile(I);
92+
System.out.println(fI.apply());
9793
test_2D_verifiy();
9894
}
9995

10096
public static void test_2D_verifiy() {
101-
double xMin=-1, xMax=1, xStep=0.001;
102-
double yMin=-1, yMax=1, yStep=0.001;
97+
double xMin=0, xMax=1, xStep=0.001;
98+
double yMin=0, yMax=1, yStep=0.001;
10399
double sum = 0.0;
104100
for(double x=xMin; x<=xMax; x+=xStep) {
105101
for(double y=yMin; y<=yMax; y+=yStep) {
@@ -138,8 +134,68 @@ public static void test_ND() {
138134
System.out.println(0.5*Math.PI*Math.PI*Math.pow(r, 4));
139135

140136
}
137+
138+
/**
139+
* I = \int_{\Omega} sin(sqrt(log(x+y+1))) dxdy
140+
* where
141+
* \Omega={ (x,y), where (x-1/2)^2 + (y-1/2)^2 <= 0.25 }
142+
*/
143+
public static void test_paper_example1() {
144+
Domain omega = new Domain2D("\\Omega", x, y)
145+
.setBound(x, 0.5-sqrt(0.25-(y-0.5)*(y-0.5)), 0.5+sqrt(0.25-(y-0.5)*(y-0.5)))
146+
.setBound(y, 0, 1)
147+
.setStepSize(0.001);
141148

142-
}
149+
Expr I = Integrate.apply( sin(sqrt(log(x+y+1)) ), omega);
150+
System.out.println(I);
143151

152+
BytecodeFunc fI = JIT.compile(I);
153+
System.out.println(fI.apply());
154+
}
155+
156+
/**
157+
* I = \int_{\Omega} sin(sqrt(log(x+y+1))) dxdy
158+
* where
159+
* \Omega= { (x,y), where
160+
* a^2 <= (x-1/2)^2 + (y-1/2)^2 <= b^2 or
161+
* c^2 <= (x-1/2)^2 + (y-1/2)^2 <= d^2
162+
* }
163+
* we choose
164+
* a=0.25, b=0.5, c=0.75, d=1.0
165+
*/
166+
public static void test_paper_example2() {
167+
168+
Expr eq = (x-0.5)*(x-0.5) + (y-0.5)*(y-0.5);
169+
Domain omega = new Domain2D("\\Omega", x, y)
170+
.setConstraint(
171+
( Ge.apply(eq, a*a) & Le.apply(eq, b*b)) |
172+
( Ge.apply(eq, c*c) & Le.apply(eq, d*d) )
173+
).setBound(x, 0, 1).setBound(y, 0, 1);
144174

175+
Expr I = Integrate.apply( sin(sqrt(log(x+y+1)) ), omega);
176+
System.out.println(I);
145177

178+
BytecodeFunc fI = JIT.compile(new Expr[]{a, b, c, d}, I);
179+
System.out.println(fI.apply(0.25, 0.5, 0.75, 1.0));
180+
181+
test_paper_example_verifiy();
182+
}
183+
184+
public static void test_paper_example_verifiy() {
185+
double xMin=0, xMax=1, xStep=0.001;
186+
double yMin=0, yMax=1, yStep=0.001;
187+
double sum = 0.0;
188+
double a=0.25, b=0.5, c=0.75, d=1.0;
189+
for(double x=xMin; x<=xMax; x+=xStep) {
190+
for(double y=yMin; y<=yMax; y+=yStep) {
191+
double disk = (x-0.5)*(x-0.5) + (y-0.5)*(y-0.5);
192+
if((a*a <= disk && disk <= b*b) ||
193+
(c*c <= disk && disk <= d*d )
194+
) {
195+
sum += Math.sin(Math.sqrt(Math.log(x+y+1)))*xStep*yStep;
196+
}
197+
}
198+
}
199+
System.out.println("verify="+sum);
200+
}
201+
}
Collapse file

‎src/symjava/relational/Ge.java‎

Copy file name to clipboardExpand all lines: src/symjava/relational/Ge.java
+1-2Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,7 @@
44
import symjava.symbolic.arity.BinaryOp;
55

66
public class Ge extends BinaryOp implements Relation {
7-
public static Expr stackTop = null;
8-
7+
98
public Ge(Expr arg1, Expr arg2) {
109
super(arg1, arg2);
1110
this.label = arg1 + " >= " + arg2;
Collapse file

‎src/symjava/symbolic/Expr.java‎

Copy file name to clipboardExpand all lines: src/symjava/symbolic/Expr.java
+12-12Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -373,18 +373,18 @@ public Expr xor(Expr other) {
373373
return Xor.simplifiedIns(this, other);
374374
}
375375

376-
/**
377-
* TODO We cannot use the comparison operator overload in java-oo for our use case
378-
* @param other
379-
* @return
380-
*/
381-
public int compareTo(Expr other) {
382-
if(Ge.stackTop == null)
383-
Ge.stackTop = Ge.apply(this, other);
384-
else
385-
Ge.stackTop = Ge.apply(Ge.stackTop, other);
386-
return -1;
387-
}
376+
// /**
377+
// * TODO We cannot use the comparison operator overload in java-oo for our use case
378+
// * @param other
379+
// * @return
380+
// */
381+
// public int compareTo(Expr other) {
382+
// if(Ge.stackTop == null)
383+
// Ge.stackTop = Ge.apply(this, other); //fix: use push()
384+
// else
385+
// Ge.stackTop = Ge.apply(Ge.stackTop, other);
386+
// return -1;
387+
// }
388388

389389

390390
/**
Collapse file

‎src/symjava/symbolic/utils/BytecodeUtils.java‎

Copy file name to clipboardExpand all lines: src/symjava/symbolic/utils/BytecodeUtils.java
+8-2Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -349,8 +349,8 @@ public static ClassGen genClass(Func fun, boolean writeClassFile, boolean static
349349
Expr x = coord[0];
350350
Expr y = coord[1];
351351
Expr xMin = INT.domain.getMinBound(x);
352-
Expr xMax = INT.domain.getMaxBound(y);
353-
//Expr yMin = INT.domain.getMinBound(x);
352+
Expr xMax = INT.domain.getMaxBound(x);
353+
//Expr yMin = INT.domain.getMinBound(y);
354354
//Expr yMax = INT.domain.getMaxBound(y);
355355
Func fxMin = new Func("integrate_bound_"+x+"Min_"+java.util.UUID.randomUUID().toString().replaceAll("-", ""), xMin);
356356
Func fxMax = new Func("integrate_bound_"+x+"Max_"+java.util.UUID.randomUUID().toString().replaceAll("-", ""), xMax);
@@ -364,6 +364,12 @@ public static ClassGen genClass(Func fun, boolean writeClassFile, boolean static
364364
//fyMin.toBytecodeFunc(true, true);
365365
//fyMax.toBytecodeFunc(true, true);
366366
//We have begin,end parameters on the top of the VM stack
367+
if(INT.domain.getStepSize(x) == null) {
368+
throw new RuntimeException("Please specify step size for "+x);
369+
}
370+
if(INT.domain.getStepSize(y) == null) {
371+
throw new RuntimeException("Please specify step size for "+y);
372+
}
367373
il.append(new PUSH(cp, INT.domain.getStepSize(y)));
368374
il.append(new PUSH(cp, fxMin.getName()));
369375
il.append(new PUSH(cp, fxMax.getName()));

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.