Skip to content

Commit cc3cc0c

Browse files
jschuellerjschueller
committed
Add RotatedCopula
1 parent 0889453 commit cc3cc0c

File tree

10 files changed

+331
-0
lines changed

10 files changed

+331
-0
lines changed

‎lib/src/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,11 @@
22

33
ot_add_current_dir_to_include_dirs ()
44

5+
ot_add_source_file (RotatedCopula.cxx)
56
ot_add_source_file (VineCopulaFactory.cxx)
67
ot_add_source_file (VineCopula.cxx)
78

9+
ot_install_header_file (RotatedCopula.hxx)
810
ot_install_header_file (VineCopulaFactory.hxx)
911
ot_install_header_file (VineCopula.hxx)
1012

‎lib/src/RotatedCopula.cxx‎

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
// -*- C++ -*-
2+
/**
3+
* @brief RotatedCopula
4+
*
5+
* Copyright 2005-2025 Airbus-EDF-IMACS-ONERA-Phimeca
6+
*
7+
* This library is free software: you can redistribute it and/or modify
8+
* it under the terms of the GNU Lesser General Public License as published by
9+
* the Free Software Foundation, either version 3 of the License, or
10+
* (at your option) any later version.
11+
*
12+
* This library is distributed in the hope that it will be useful,
13+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
14+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15+
* GNU Lesser General Public License for more details.
16+
*
17+
* You should have received a copy of the GNU Lesser General Public License
18+
* along with this library. If not, see <http://www.gnu.org/licenses/>.
19+
*
20+
*/
21+
#include "otvine/RotatedCopula.hxx"
22+
#include <openturns/PersistentObjectFactory.hxx>
23+
24+
#include <vinecopulib.hpp>
25+
26+
using namespace OT;
27+
28+
namespace OTVINE
29+
{
30+
31+
CLASSNAMEINIT(RotatedCopula);
32+
33+
static Factory<RotatedCopula> Factory_RotatedCopula;
34+
35+
36+
/* Default constructor */
37+
RotatedCopula::RotatedCopula()
38+
: DistributionImplementation()
39+
{
40+
isCopula_ = true;
41+
setName("RotatedCopula");
42+
}
43+
44+
RotatedCopula::RotatedCopula(const OT::Distribution & copula, const Rotation rotation)
45+
: DistributionImplementation()
46+
, copula_(copula)
47+
, rotation_(rotation)
48+
{
49+
isCopula_ = true;
50+
setName("RotatedCopula");
51+
setDimension(copula.getDimension());
52+
setDescription(copula.getDescription());
53+
computeRange();
54+
if (copula.getDimension() != 2)
55+
throw NotYetImplementedException(HERE) << "RotatedCopula only supports dimension 2";
56+
}
57+
58+
/* Virtual constructor method */
59+
RotatedCopula * RotatedCopula::clone() const
60+
{
61+
return new RotatedCopula(*this);
62+
}
63+
64+
Point RotatedCopula::rotate(const Point & x, const Bool inverse) const
65+
{
66+
if (x.getDimension() != 2)
67+
throw NotYetImplementedException(HERE) << "RotatedCopula only supports dimension 2";
68+
const Scalar u1 = x[0];
69+
const Scalar u2 = x[1];
70+
Rotation rotation = static_cast<Rotation>(rotation_);
71+
if (inverse)
72+
{
73+
switch (rotation)
74+
{
75+
case (ROTATION_90):
76+
rotation = ROTATION_270;
77+
break;
78+
case (ROTATION_270):
79+
rotation = ROTATION_90;
80+
break;
81+
default:
82+
break;
83+
}
84+
}
85+
86+
switch (rotation)
87+
{
88+
case (ROTATION_0):
89+
return {u1, u2};
90+
case (ROTATION_90):
91+
return {u2, 1.0 - u1};
92+
case (ROTATION_180):
93+
return {1.0 - u2, 1.0 - u1};
94+
case (ROTATION_270):
95+
return {1.0 - u2, u1};
96+
default:
97+
throw NotYetImplementedException(HERE) << "RotatedCopula invalid rotation";
98+
}
99+
}
100+
101+
/* Get the PDF of the distribution */
102+
Scalar RotatedCopula::computePDF(const Point & point) const
103+
{
104+
return copula_.computePDF(rotate(point));
105+
}
106+
107+
/* Get the CDF of the distribution */
108+
Scalar RotatedCopula::computeCDF(const Point & point) const
109+
{
110+
// https://quantitative-thinking.com/copula_theory/cop_modifications/clockwise_rot.html
111+
const Scalar cdf0 = copula_.computeCDF(rotate(point));
112+
switch (rotation_)
113+
{
114+
case (ROTATION_0):
115+
return cdf0;
116+
case (ROTATION_90):
117+
return copula_.computeCDF(Point({point[1], 1.0})) - cdf0;
118+
case (ROTATION_180):
119+
return 1.0 - copula_.computeCDF(Point({1.0, 1.0 - point[1]})) - copula_.computeCDF(Point({1.0 - point[0], 1.0})) + 2.0 * cdf0;
120+
case (ROTATION_270):
121+
return copula_.computeCDF(Point({1.0, point[0]})) - cdf0;
122+
default:
123+
throw NotYetImplementedException(HERE) << "RotatedCopula invalid rotation";
124+
}
125+
}
126+
127+
/* Get one realization of the distribution */
128+
Point RotatedCopula::getRealization() const
129+
{
130+
const Point x = copula_.getRealization();
131+
return rotate(x, true);
132+
}
133+
134+
/* Parameters value accessor */
135+
Point RotatedCopula::getParameter() const
136+
{
137+
return copula_.getParameter();
138+
}
139+
140+
void RotatedCopula::setParameter(const Point & parameter)
141+
{
142+
copula_.setParameter(parameter);
143+
}
144+
145+
/* Parameters description accessor */
146+
Description RotatedCopula::getParameterDescription() const
147+
{
148+
return copula_.getParameterDescription();
149+
}
150+
151+
/* Comparison operator */
152+
Bool RotatedCopula::operator ==(const RotatedCopula & other) const
153+
{
154+
if (this == &other) return true;
155+
return (copula_ == other.copula_) && (rotation_ == other.rotation_);
156+
}
157+
158+
Bool RotatedCopula::equals(const DistributionImplementation & other) const
159+
{
160+
const RotatedCopula* p_other = dynamic_cast<const RotatedCopula*>(&other);
161+
return p_other && (*this == *p_other);
162+
}
163+
164+
/* String converter */
165+
String RotatedCopula::__repr__() const
166+
{
167+
OSS oss;
168+
oss << "class=" << RotatedCopula::GetClassName() << " copula=" << copula_ << " rotation=" << rotation_;
169+
return oss;
170+
}
171+
172+
/* Method save() stores the object through the StorageManager */
173+
void RotatedCopula::save(Advocate & adv) const
174+
{
175+
DistributionImplementation::save(adv);
176+
adv.saveAttribute("copula_", copula_);
177+
adv.saveAttribute("rotation_", rotation_);
178+
}
179+
180+
/* Method load() reloads the object from the StorageManager */
181+
void RotatedCopula::load(Advocate & adv)
182+
{
183+
DistributionImplementation::load(adv);
184+
adv.loadAttribute("copula_", copula_);
185+
adv.loadAttribute("rotation_", rotation_);
186+
computeRange();
187+
}
188+
189+
190+
} /* namespace OTVINE */

‎lib/src/otvine/RotatedCopula.hxx‎

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
// -*- C++ -*-
2+
/**
3+
* @brief RotatedCopula
4+
*
5+
* Copyright 2005-2025 Airbus-EDF-IMACS-ONERA-Phimeca
6+
*
7+
* This library is free software: you can redistribute it and/or modify
8+
* it under the terms of the GNU Lesser General Public License as published by
9+
* the Free Software Foundation, either version 3 of the License, or
10+
* (at your option) any later version.
11+
*
12+
* This library is distributed in the hope that it will be useful,
13+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
14+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15+
* GNU Lesser General Public License for more details.
16+
*
17+
* You should have received a copy of the GNU Lesser General Public License
18+
* along with this library. If not, see <http://www.gnu.org/licenses/>.
19+
*
20+
*/
21+
#ifndef OTVINE_ROTATEDCOPULA_HXX
22+
#define OTVINE_ROTATEDCOPULA_HXX
23+
24+
#include <openturns/DistributionImplementation.hxx>
25+
#include <openturns/Distribution.hxx>
26+
27+
#include "otvine/otvineprivate.hxx"
28+
29+
namespace vinecopulib {
30+
class Vinecop;
31+
}
32+
33+
namespace OTVINE
34+
{
35+
36+
/**
37+
* @class RotatedCopula
38+
*
39+
* RotatedCopula is some vinecopulafactory type to illustrate how to add some classes in OpenTURNS
40+
*/
41+
class OTVINE_API RotatedCopula
42+
: public OT::DistributionImplementation
43+
{
44+
CLASSNAME
45+
46+
public:
47+
/** Default constructor */
48+
RotatedCopula();
49+
50+
enum Rotation {ROTATION_0, ROTATION_90, ROTATION_180, ROTATION_270};
51+
52+
explicit RotatedCopula(const OT::Distribution & copula, const Rotation = ROTATION_90);
53+
54+
/** Virtual constructor method */
55+
RotatedCopula * clone() const override;
56+
57+
OT::Bool operator ==(const RotatedCopula & other) const;
58+
protected:
59+
OT::Bool equals(const DistributionImplementation & other) const override;
60+
public:
61+
62+
/** String converter */
63+
OT::String __repr__() const override;
64+
65+
/** Get one realization of the distribution */
66+
OT::Point getRealization() const override;
67+
68+
/** Get the PDF of the distribution */
69+
using OT::DistributionImplementation::computePDF;
70+
OT::Scalar computePDF(const OT::Point & point) const override;
71+
72+
/** Get the CDF of the distribution */
73+
using OT::DistributionImplementation::computeCDF;
74+
OT::Scalar computeCDF(const OT::Point & point) const override;
75+
76+
/** Parameters value accessors */
77+
void setParameter(const OT::Point & parameter) override;
78+
OT::Point getParameter() const override;
79+
80+
/** Parameters description accessor */
81+
OT::Description getParameterDescription() const override;
82+
83+
/** Method save() stores the object through the StorageManager */
84+
void save(OT::Advocate & adv) const override;
85+
86+
/** Method load() reloads the object from the StorageManager */
87+
void load(OT::Advocate & adv) override;
88+
89+
private:
90+
OT::Point rotate(const OT::Point & x, const OT::Bool inverse = false) const;
91+
92+
OT::Distribution copula_;
93+
OT::UnsignedInteger rotation_ = ROTATION_0;
94+
95+
}; /* class RotatedCopula */
96+
97+
} /* namespace OTVINE */
98+
99+
#endif /* OTVINE_ROTATEDCOPULA_HXX */

‎python/doc/user_manual/user_manual.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,5 @@ API Reference
88
:template: class.rst_t
99

1010
VineCopulaFactory
11+
VineCopula
12+
RotatedCopula

‎python/src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ endmacro ()
4747

4848

4949
ot_add_python_module( ${PACKAGE_NAME} ${PACKAGE_NAME}_module.i
50+
RotatedCopula.i RotatedCopula_doc.i
5051
VineCopula.i VineCopula_doc.i
5152
VineCopulaFactory.i VineCopulaFactory_doc.i
5253
)

‎python/src/RotatedCopula.i‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
// SWIG file RotatedCopula.i
2+
3+
%{
4+
#include "otvine/RotatedCopula.hxx"
5+
%}
6+
7+
%include RotatedCopula_doc.i
8+
9+
%copyctor OTVINE::RotatedCopula;
10+
11+
%include otvine/RotatedCopula.hxx

‎python/src/RotatedCopula_doc.i‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
%feature("docstring") OTVINE::RotatedCopula
2+
"Rotated copula.
3+
4+
Examples
5+
--------
6+
>>> import otvine
7+
>>> import openturns as ot
8+
>>> copula = ot.GumbelCopula()
9+
>>> rotated = otvine.RotatedCopula(copula)
10+
"

‎python/src/otvine_module.i‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818

1919
// The new classes
2020
%include otvine/otvineprivate.hxx
21+
%include RotatedCopula.i
2122
%include VineCopula.i
2223
%include VineCopulaFactory.i
2324

‎python/test/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ macro (ot_pyinstallcheck_test TESTNAME)
3737
endmacro ()
3838

3939

40+
ot_pyinstallcheck_test (RotatedCopula_std IGNOREOUT)
4041
ot_pyinstallcheck_test (VineCopulaFactory_std IGNOREOUT)
4142
ot_pyinstallcheck_test (docstring IGNOREOUT)
4243

‎python/test/t_RotatedCopula_std.py‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
#!/usr/bin/env python
2+
3+
import openturns as ot
4+
import openturns.testing as ott
5+
import otvine
6+
7+
copula = ot.GumbelCopula()
8+
rotated = otvine.RotatedCopula(copula)
9+
10+
print(ot.RandomGenerator.Generate())
11+
12+
ot.Log.Show(ot.Log.TRACE)
13+
validation = ott.DistributionValidation(rotated)
14+
validation.run()

0 commit comments

Comments
 (0)