From 19929cfdd66500e7d2fa4655550b8004abe42174 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 15 Oct 2024 22:22:23 +0200 Subject: [PATCH 1/6] add skeleton for compiled support package based on pybind11 and ducc0 --- .gitmodules | 3 + external/ducc0 | 1 + native_support/cmdr4_support.cc | 28 ++++++ native_support/pyproject.toml | 24 +++++ native_support/python/utils_pymod.cc | 70 ++++++++++++++ native_support/setup.py | 132 +++++++++++++++++++++++++++ 6 files changed, 258 insertions(+) create mode 100644 .gitmodules create mode 160000 external/ducc0 create mode 100644 native_support/cmdr4_support.cc create mode 100644 native_support/pyproject.toml create mode 100644 native_support/python/utils_pymod.cc create mode 100644 native_support/setup.py diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..b6b0201 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "ducc0"] + path = external/ducc0 + url = https://gitlab.mpcdf.mpg.de/mtr/ducc.git diff --git a/external/ducc0 b/external/ducc0 new file mode 160000 index 0000000..e0ef9c2 --- /dev/null +++ b/external/ducc0 @@ -0,0 +1 @@ +Subproject commit e0ef9c25fc74b78392359745cbb39b0b445cb7ea diff --git a/native_support/cmdr4_support.cc b/native_support/cmdr4_support.cc new file mode 100644 index 0000000..3558857 --- /dev/null +++ b/native_support/cmdr4_support.cc @@ -0,0 +1,28 @@ +#include "ducc0/infra/string_utils.cc" +#include "ducc0/infra/threading.cc" +#include "ducc0/infra/mav.cc" +//#include "ducc0/math/pointing.cc" +//#include "ducc0/math/geom_utils.cc" +//#include "ducc0/math/space_filling.cc" +//#include "ducc0/math/gl_integrator.cc" +//#include "ducc0/math/gridding_kernel.cc" +//#include "ducc0/math/wigner3j.cc" +//#include "ducc0/sht/sht.cc" +//#include "ducc0/healpix/healpix_tables.cc" +//#include "ducc0/healpix/healpix_base.cc" + +#include +#include "python/utils_pymod.cc" + +using namespace cmdr4; + +PYBIND11_MODULE(PKGNAME, m) + { +#define CMDR4_XSTRINGIFY(s) CMDR4_STRINGIFY(s) +#define CMDR4_STRINGIFY(s) #s + m.attr("__version__") = CMDR4_XSTRINGIFY(PKGVERSION); +#undef CMDR4_STRINGIFY +#undef CMDR4_XSTRINGIFY + + add_utils(m); + } diff --git a/native_support/pyproject.toml b/native_support/pyproject.toml new file mode 100644 index 0000000..caa5bf2 --- /dev/null +++ b/native_support/pyproject.toml @@ -0,0 +1,24 @@ +[project] +name = "cmdr4_support" +description = "Compiled support code for Commander4" +authors = [ + {name = "Martin Reinecke", email = "martin@mpa-garching.mpg.de"}, +] +readme = "README.md" +license = {file = "LICENSE"} +classifiers = [ + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Mathematics", + "Topic :: Scientific/Engineering :: Physics", + "License :: OSI Approved :: GNU General Public License v2 or later (GPLv2+)", + "Operating System :: OS Independent", + "Programming Language :: C++", + "Programming Language :: Python", +] +requires-python = ">=3.8" +dependencies = ["numpy>=1.17.0"] +dynamic = ["version"] + +[build-system] +requires = ["setuptools >= 40.6.0", "pybind11 >= 2.6.0", "numpy >= 1.17.0"] +build-backend = "setuptools.build_meta" diff --git a/native_support/python/utils_pymod.cc b/native_support/python/utils_pymod.cc new file mode 100644 index 0000000..e550736 --- /dev/null +++ b/native_support/python/utils_pymod.cc @@ -0,0 +1,70 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "ducc0/infra/mav.h" +#include "ducc0/infra/misc_utils.h" +#include "ducc0/math/constants.h" +#include "ducc0/bindings/pybind_utils.h" + +namespace cmdr4 { + +namespace detail_pymodule_utils { + +using namespace std; +using namespace ducc0; +namespace py = pybind11; +auto None = py::none(); + + +py::array Py_amplitude_sampling_per_pix_helper ( + const py::array &map_sky_, + const py::array &map_rms_, + const py::array &M_, + py::object &comp_maps__, + size_t nthreads) + { + const auto map_sky = to_cmav(map_sky_); + size_t nband = map_sky.shape(0); + size_t npix = map_sky.shape(1); + const auto map_rms = to_cmav(map_rms_); + MR_assert(map_rms.shape()==map_sky.shape(), "map shape mismatch"); + const auto M = to_cmav(M_); + size_t ncomp = M.shape(1); + MR_assert(M.shape(0)==nband, "M shape mismatch"); + auto comp_maps_ = get_optional_Pyarr(comp_maps__, {ncomp, npix}); + auto comp_maps = to_vmav(comp_maps_); + { + py::gil_scoped_release release; + +/* do computations here */ + + } + return comp_maps_; + } + +void add_utils(py::module_ &msup) + { + using namespace pybind11::literals; + auto m = msup.def_submodule("utils"); +// m.doc() = utils_DS; + + m.def("amplitude_sampling_per_pix_helper", + Py_amplitude_sampling_per_pix_helper, +// Py_amplitude_sampling_per_pix_helper_DS, + "map_sky"_a, + "map_rms"_a, + "M"_a, + "comp_maps"_a=None, + "nthreads"_a=1); + } + +} + +using detail_pymodule_utils::add_utils; + +} diff --git a/native_support/setup.py b/native_support/setup.py new file mode 100644 index 0000000..b54effd --- /dev/null +++ b/native_support/setup.py @@ -0,0 +1,132 @@ +import sys +import os.path +import itertools +from glob import iglob +import os + +from setuptools import setup, Extension +import pybind11 + +pkgname = 'cmdr4_support' +version = '0.0.1' + +tmp = os.getenv("CMDR4_CFLAGS", "").split(" ") +user_cflags = [x for x in tmp if x != ""] +tmp = os.getenv("CMDR4_LFLAGS", "").split(" ") +user_lflags = [x for x in tmp if x != ""] +tmp = os.getenv("CMDR4_FLAGS", "").split(" ") +tmp = [x for x in tmp if x != ""] +user_cflags += tmp +user_lflags += tmp + +compilation_strategy = os.getenv('CMDR4_OPTIMIZATION', 'native-strip') +if compilation_strategy not in ['none', 'none-debug', 'none-strip', 'portable', 'portable-debug', 'portable-strip', 'native', 'native-debug', 'native-strip']: + raise RuntimeError('unknown compilation strategy') +do_debug = compilation_strategy in ['none-debug', 'portable-debug', 'native-debug'] +do_strip = compilation_strategy in ['none-strip', 'portable-strip', 'native-strip'] +do_optimize = compilation_strategy not in ['none', 'none-debug', 'none-strip'] +do_native = compilation_strategy in ['native', 'native-debug', 'native-strip'] + +def _print_env(): + import platform + print("") + print("Build environment:") + print("Platform: ", platform.platform()) + print("Machine: ", platform.machine()) + print("System: ", platform.system()) + print("Architecture: ", platform.architecture()) + print("") + +def _get_files_by_suffix(directory, suffix): + path = directory + iterable_sources = (iglob(os.path.join(root, '*.'+suffix)) + for root, dirs, files in os.walk(path)) + return list(itertools.chain.from_iterable(iterable_sources)) + + +include_dirs = ['.', '../external/ducc0/src/', + pybind11.get_include(True), + pybind11.get_include(False)] + +extra_compile_args = ['-std=c++17', '-fvisibility=hidden'] + +if do_debug: + extra_compile_args += ['-g'] +else: + extra_compile_args += ['-g0'] + +if do_optimize: + extra_compile_args += ['-ffast-math', +# '-ffp-contract=fast', +# '-fno-math-errno', +# '-freciprocal-math', +# '-fassociative-math', +# '-fno-signed-zeros', +# '-fno-trapping-math', +# '-fno-signaling-nans', +# '-fcx-limited-range', +# '-fexcess-precision=fast', +# '-ffinite-math-only', + '-O3'] +else: + extra_compile_args += ['-O0'] + +if do_native: + extra_compile_args += ['-march=native'] + +python_module_link_args = [] + +define_macros = [("PKGNAME", pkgname), + ("PKGVERSION", version)] + +if sys.platform == 'darwin': + extra_compile_args += ['-mmacosx-version-min=10.14', '-pthread'] + python_module_link_args += ['-mmacosx-version-min=10.14', '-pthread'] +elif sys.platform == 'win32': + extra_compile_args = ['/EHsc', '/std:c++17'] + if do_optimize: + extra_compile_args += ['/Ox'] +else: + if do_optimize: + extra_compile_args += ['-Wfatal-errors', + '-Wfloat-conversion', + '-W', + '-Wall', + '-Wstrict-aliasing', + '-Wwrite-strings', + '-Wredundant-decls', + '-Woverloaded-virtual', + '-Wcast-qual', + '-Wcast-align', + '-Wpointer-arith', + '-Wnon-virtual-dtor', + '-Wzero-as-null-pointer-constant'] + extra_compile_args += ['-pthread'] + python_module_link_args += ['-Wl,-rpath,$ORIGIN', '-pthread'] + if do_native: + python_module_link_args += ['-march=native'] + if do_strip: + python_module_link_args += ['-s'] + +extra_compile_args += user_cflags +python_module_link_args += user_lflags + +depfiles = (_get_files_by_suffix('.', 'h') + + _get_files_by_suffix('.', 'cc') + + ['setup.py']) + +extensions = [Extension(pkgname, + language='c++', + sources=['cmdr4_support.cc'], + depends=depfiles, + include_dirs=include_dirs, + define_macros=define_macros, + extra_compile_args=extra_compile_args, + extra_link_args=python_module_link_args)] + +_print_env() + +setup(name=pkgname, + version=version, + ext_modules = extensions + ) From 3d386c780313ba2917aba19363e060b879b5b67d Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Thu, 17 Oct 2024 11:17:26 +0200 Subject: [PATCH 2/6] implement native amplitude sampling helper --- native_support/python/utils_pymod.cc | 205 +++++++++++++++++++++++++++ src/python/compsep_loop.py | 12 +- src/python/tod_loop.py | 4 +- 3 files changed, 218 insertions(+), 3 deletions(-) diff --git a/native_support/python/utils_pymod.cc b/native_support/python/utils_pymod.cc index e550736..22bc43b 100644 --- a/native_support/python/utils_pymod.cc +++ b/native_support/python/utils_pymod.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -20,11 +21,169 @@ using namespace ducc0; namespace py = pybind11; auto None = py::none(); +#if 1 + +struct SingularError {}; + +//! Find pivot element +/*! +* The function pivot finds the largest element for a pivot in "jcol" +* of Matrix "a", performs interchanges of the appropriate +* rows in "a", and also interchanges the corresponding elements in +* the order vector. +* +* +* \param a - n by n Matrix of coefficients +* \param order - integer vector to hold row ordering +* \param jcol - column of "a" being searched for pivot element +* +*/ +int pivot(vmav &a, vmav &order, int jcol) + { + constexpr double TINY=1e-20; + int n = a.shape(0); + + /* + * Find biggest element on or below diagonal. + * This will be the pivot row. + */ + + int ipvt = jcol; + double big = std::abs(a(ipvt,ipvt)); + for (int i = ipvt+1; ibig) + { + big = anext; + ipvt = i; + } + } + if (std::abs(big)TINY); // otherwise Matrix is singular + + /* Interchange pivot row (ipvt) with current row (jcol). */ + + if (ipvt==jcol) return 0; + //a.swaprows(jcol,ipvt); + for (int i=0; i &a, vmav &order) + { + int flag = 1; /* changes sign with each row interchange */ + + int n = a.shape(0); + MR_assert(a.shape(1)==size_t(n)); + MR_assert(order.shape(0)==size_t(n)); + + /* establish initial ordering in order vector */ + for (int i=0; i &a, const cmav &b, vmav &x, const cmav &order) + { + int n = a.shape(0); + + /* rearrange the elements of the b vector. x is used to hold them. */ + + for (int i=0; i=0; i--) + { + double sum = 0.0; + for (int j=i+1; j(M_); size_t ncomp = M.shape(1); MR_assert(M.shape(0)==nband, "M shape mismatch"); + const auto random = to_cmav(random_); + MR_assert(random.shape(0)==npix, "random numbers shape mismatch"); + MR_assert(random.shape(1)==nband, "random numbers shape mismatch"); auto comp_maps_ = get_optional_Pyarr(comp_maps__, {ncomp, npix}); auto comp_maps = to_vmav(comp_maps_); { py::gil_scoped_release release; /* do computations here */ + execStatic(npix, nthreads, 0,[&](auto &sched) { + + vmav x({ncomp}), xmap({nband}), tmap({nband}), comp({ncomp}); + vmav order({ncomp}); + vmav A({ncomp,ncomp}); + while (auto rng=sched.getNext()) + for (size_t it=rng.lo; it Date: Thu, 17 Oct 2024 12:16:31 +0200 Subject: [PATCH 3/6] cleanup --- native_support/python/utils_pymod.cc | 34 ++++++++-------------------- 1 file changed, 10 insertions(+), 24 deletions(-) diff --git a/native_support/python/utils_pymod.cc b/native_support/python/utils_pymod.cc index 22bc43b..768f37a 100644 --- a/native_support/python/utils_pymod.cc +++ b/native_support/python/utils_pymod.cc @@ -21,8 +21,6 @@ using namespace ducc0; namespace py = pybind11; auto None = py::none(); -#if 1 - struct SingularError {}; //! Find pivot element @@ -61,10 +59,8 @@ int pivot(vmav &a, vmav &order, int jcol) } if (std::abs(big)TINY); // otherwise Matrix is singular - + /* Interchange pivot row (ipvt) with current row (jcol). */ - if (ipvt==jcol) return 0; //a.swaprows(jcol,ipvt); for (int i=0; i &a, const cmav &b, vmav &x, int n = a.shape(0); /* rearrange the elements of the b vector. x is used to hold them. */ - for (int i=0; i &a, const cmav &b, vmav &x, x(i) -= sum; } } -#endif py::array Py_amplitude_sampling_per_pix_helper ( const py::array &map_sky_, @@ -203,9 +193,8 @@ py::array Py_amplitude_sampling_per_pix_helper ( { py::gil_scoped_release release; -/* do computations here */ - execStatic(npix, nthreads, 0,[&](auto &sched) { - + execStatic(npix, nthreads, 0,[&](auto &sched) + { vmav x({ncomp}), xmap({nband}), tmap({nband}), comp({ncomp}); vmav order({ncomp}); vmav A({ncomp,ncomp}); @@ -221,32 +210,29 @@ py::array Py_amplitude_sampling_per_pix_helper ( { x(i) = 0; for (size_t j=0; j Date: Thu, 17 Oct 2024 12:31:51 +0200 Subject: [PATCH 4/6] add docstring --- native_support/python/utils_pymod.cc | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/native_support/python/utils_pymod.cc b/native_support/python/utils_pymod.cc index 768f37a..36ec846 100644 --- a/native_support/python/utils_pymod.cc +++ b/native_support/python/utils_pymod.cc @@ -237,6 +237,30 @@ py::array Py_amplitude_sampling_per_pix_helper ( return comp_maps_; } +constexpr const char *Py_amplitude_sampling_per_pix_helper_DS = R"""( +Helper function for amplitude_ampling_per_pix. + +Parameters +---------- +map_sky: numpy.ndarray((nband, npix), dtype=np.float64) + the sky maps for all bands +map_rms: numpy.ndarray((nband, npix), dtype=np.float64) + the RMS maps for all bands +M: numpy.ndarray((nband, ncomp), dtype=np.float64) + the mixing matrix +random: numpy.ndarray((npix, nband), dtype=np.float64) + normal-distributed random numbers + NOTE: the transposed shape is only for performance reasons +comp_maps: numpy.ndarray((npix, ncomp), dtype=np.float64) or None + If provided, the results will be stored in this array + +Returns +------- +numpy.ndarray((npix, ncomp), dtype=np.float64) + The component maps + If comp_maps was provided, this array is identical to comp_maps +)"""; + void add_utils(py::module_ &msup) { using namespace pybind11::literals; @@ -245,7 +269,7 @@ void add_utils(py::module_ &msup) m.def("amplitude_sampling_per_pix_helper", Py_amplitude_sampling_per_pix_helper, -// Py_amplitude_sampling_per_pix_helper_DS, + Py_amplitude_sampling_per_pix_helper_DS, "map_sky"_a, "map_rms"_a, "M"_a, From a11f8fc0fee94b0ccbb6482ecab0b6635747aa8f Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Fri, 18 Oct 2024 14:33:11 +0200 Subject: [PATCH 5/6] cosmetics --- README.md | 17 ++++++++++- cmdr4_support/cmdr4_support.cc | 15 ++++++++++ .../pyproject.toml | 0 .../python/utils_pymod.cc | 1 + {native_support => cmdr4_support}/setup.py | 0 native_support/cmdr4_support.cc | 28 ------------------- src/python/compsep_loop.py | 4 ++- 7 files changed, 35 insertions(+), 30 deletions(-) create mode 100644 cmdr4_support/cmdr4_support.cc rename {native_support => cmdr4_support}/pyproject.toml (100%) rename {native_support => cmdr4_support}/python/utils_pymod.cc (99%) rename {native_support => cmdr4_support}/setup.py (100%) delete mode 100644 native_support/cmdr4_support.cc diff --git a/README.md b/README.md index 53364cc..de5c225 100644 --- a/README.md +++ b/README.md @@ -8,4 +8,19 @@ mpirun -n 15 python -u src/python/commander4.py -p params/param_default.yml (The `-u` makes stdout unbuffered, which I have found to be necessary in order to make MPI programs print properly to the terminal). ## Compiling Ctypes libraries -The code depends on C++ Ctypes libraries which are located in the `src/cpp/` directory. There should be a file named `src/cpp/compilation.sh` which contains the compilation procedure for all C++ files. \ No newline at end of file +The code depends on C++ Ctypes libraries which are located in the `src/cpp/` directory. There should be a file named `src/cpp/compilation.sh` which contains the compilation procedure for all C++ files. + +## Initializing the submodules +Commander4 pulls in the ducc0 sources to make developing of C++ helper functions easier. + +To install this submodule directly when cloning the Commander4 repository you can do +``` +git clone --recurse-submodules +``` + +If you have already cloned Commander 4, the easiest way is to go to the Commender4 directory +and then do +``` +git submodule init +git submodule update +``` diff --git a/cmdr4_support/cmdr4_support.cc b/cmdr4_support/cmdr4_support.cc new file mode 100644 index 0000000..fb88dd0 --- /dev/null +++ b/cmdr4_support/cmdr4_support.cc @@ -0,0 +1,15 @@ +#include +#include "python/utils_pymod.cc" + +using namespace cmdr4; + +PYBIND11_MODULE(PKGNAME, m) + { +#define CMDR4_XSTRINGIFY(s) CMDR4_STRINGIFY(s) +#define CMDR4_STRINGIFY(s) #s + m.attr("__version__") = CMDR4_XSTRINGIFY(PKGVERSION); +#undef CMDR4_STRINGIFY +#undef CMDR4_XSTRINGIFY + + add_utils(m); + } diff --git a/native_support/pyproject.toml b/cmdr4_support/pyproject.toml similarity index 100% rename from native_support/pyproject.toml rename to cmdr4_support/pyproject.toml diff --git a/native_support/python/utils_pymod.cc b/cmdr4_support/python/utils_pymod.cc similarity index 99% rename from native_support/python/utils_pymod.cc rename to cmdr4_support/python/utils_pymod.cc index 36ec846..58dd887 100644 --- a/native_support/python/utils_pymod.cc +++ b/cmdr4_support/python/utils_pymod.cc @@ -21,6 +21,7 @@ using namespace ducc0; namespace py = pybind11; auto None = py::none(); +// exception type for signalling a (nearly) singular matrix) struct SingularError {}; //! Find pivot element diff --git a/native_support/setup.py b/cmdr4_support/setup.py similarity index 100% rename from native_support/setup.py rename to cmdr4_support/setup.py diff --git a/native_support/cmdr4_support.cc b/native_support/cmdr4_support.cc deleted file mode 100644 index 3558857..0000000 --- a/native_support/cmdr4_support.cc +++ /dev/null @@ -1,28 +0,0 @@ -#include "ducc0/infra/string_utils.cc" -#include "ducc0/infra/threading.cc" -#include "ducc0/infra/mav.cc" -//#include "ducc0/math/pointing.cc" -//#include "ducc0/math/geom_utils.cc" -//#include "ducc0/math/space_filling.cc" -//#include "ducc0/math/gl_integrator.cc" -//#include "ducc0/math/gridding_kernel.cc" -//#include "ducc0/math/wigner3j.cc" -//#include "ducc0/sht/sht.cc" -//#include "ducc0/healpix/healpix_tables.cc" -//#include "ducc0/healpix/healpix_base.cc" - -#include -#include "python/utils_pymod.cc" - -using namespace cmdr4; - -PYBIND11_MODULE(PKGNAME, m) - { -#define CMDR4_XSTRINGIFY(s) CMDR4_STRINGIFY(s) -#define CMDR4_STRINGIFY(s) #s - m.attr("__version__") = CMDR4_XSTRINGIFY(PKGVERSION); -#undef CMDR4_STRINGIFY -#undef CMDR4_XSTRINGIFY - - add_utils(m); - } diff --git a/src/python/compsep_loop.py b/src/python/compsep_loop.py index 8696e44..e3fec5b 100644 --- a/src/python/compsep_loop.py +++ b/src/python/compsep_loop.py @@ -106,9 +106,11 @@ def amplitude_sampling_per_pix(map_sky: np.array, map_rms: np.array, freqs: np.a M[:,0] = CMB().get_sed(freqs) M[:,1] = ThermalDust().get_sed(freqs) M[:,2] = Synchrotron().get_sed(freqs) - rand = np.random.randn(npix,nband) from time import time t0 = time() + rand = np.random.randn(npix,nband) + print(f"time for random numbers: {time()-t0}s.") + t0 = time() for i in range(npix): xmap = 1/map_rms[:,i] x = M.T.dot((xmap**2*map_sky[:,i])) From b0d6e2d1fd1378be6cfa3a9a023672e56e8a52b1 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Fri, 18 Oct 2024 14:55:56 +0200 Subject: [PATCH 6/6] simplify --- cmdr4_support/cmdr4_support.cc | 2 +- cmdr4_support/{python => }/utils_pymod.cc | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) rename cmdr4_support/{python => }/utils_pymod.cc (98%) diff --git a/cmdr4_support/cmdr4_support.cc b/cmdr4_support/cmdr4_support.cc index fb88dd0..3b4248d 100644 --- a/cmdr4_support/cmdr4_support.cc +++ b/cmdr4_support/cmdr4_support.cc @@ -1,5 +1,5 @@ #include -#include "python/utils_pymod.cc" +#include "utils_pymod.cc" using namespace cmdr4; diff --git a/cmdr4_support/python/utils_pymod.cc b/cmdr4_support/utils_pymod.cc similarity index 98% rename from cmdr4_support/python/utils_pymod.cc rename to cmdr4_support/utils_pymod.cc index 58dd887..0cb546b 100644 --- a/cmdr4_support/python/utils_pymod.cc +++ b/cmdr4_support/utils_pymod.cc @@ -24,6 +24,8 @@ auto None = py::none(); // exception type for signalling a (nearly) singular matrix) struct SingularError {}; +// LU decomposition code taken from https://www.johnloomis.org/ece538/notes/2008/Matrix/ludcmp.html + //! Find pivot element /*! * The function pivot finds the largest element for a pivot in "jcol" @@ -93,7 +95,6 @@ int ludcmp(vmav &a, vmav &order) * The general plan is to compute a column of L's, then * call pivot to interchange rows, and then compute * a row of U's. */ - int nm1 = n - 1; for (int j=1; j